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Abstract 



This paper introduces a framework for simulating finite dimensional representations of (jump) diffusion sample 
paths over finite intervals, without discretisation error (exactly), in such a way that the sample path can be re- 
stored at any finite collection of time points. Within this framework we extend existing exact algorithms and 
introduce novel adaptive approaches. We consider an application of the methodology developed within this pa- 
per which allows the simulation of upper and lower bounding processes which almost surely constrain (jump) 
difi^usion sample paths to any specified tolerance. 
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1 Introduction 

A jump diffusion V : R — > R is a Markov process. We consider jump diffusions defined as the solution to a 
stochastic differential equation (SDE) of the form (denoting f- := lims|, s). 



where yS : R — > R and cr : R — > R+ denote the (instantaneous) drift and diffusion coefficients respectively, W, is a 
standard Brownian Motion and y/'" denotes a compound Poisson process, jf'^ is parameterised with (finite) jump 
intensity /I : R ^ R+, jump size coefficient v : R — > R and as jumps distributed with density fy. All coefficients 
are themselves (typically) dependent on V,. Regularity conditions are assumed to hold to ensure the existence of a 
unique non-explosive weak solution (see for instance Il28l[30l ). A discussion of conditions sufficient to allow the 
application of the methodology presented within this paper is given in Section|2] 

Diffusions and jump diffusions are widely used across a number of application areas. An extensive literature 
exists in economics and finance, spanning from the Black-Scholes model |T0"25^ '26| to the present fT5"4l. Other 
applications can be easily found within the physical |f29l and life sciences |18, 19). Motivated by such applica- 
tions, there is considerable interest in finding solutions to SDEs. In particular, we would like to be able to draw 
sample paths from the measure induced by ([T]), denoted T''. However, this is complicated as diffiision sample 
paths are infinite dimensional random variables (so at best we can hope to sample some finite dimensional subset 
of the sample path) and there doesn't typically exist an analytically tractable form for the transition density. 

Diffusion sample paths can be simulated approximately at a finite collection of time points by discretisation 
Il24ll30l . noting that as Brownian motion has a Gaussian transition density then over short intervals the transition 
density of ([T|l can be approximated by one with fixed coefficients (by a continuity argument). This can be achieved 
by breaking the interval the sample path is to be simulated over into a fine mesh (for instance, of size Af), then 
iteratively (at each mesh point) fixing the coefficients and then simulating the sample path to the next mesh point. 
For instance, in an Euler discretisation 12 111 of ([l]), the sample path is propagated between mesh points as follows. 



It is hoped the simulated sample path (generated approximately at a finite collection of mesh points) can be used 
as a proxy for an entire sample path drawn exactly from T''. More complex discretisation schemes exist (for 



dV, = j6(y,0 df + cr(y,_) AW, + dJ, 



Vo = ve R, te [0,T], 



(1) 




w.p. exp{-A(V,)At} , 
w.p. 1 - exp{-/l(y,)Af} . 



(2) 
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instance, by exploiting Ito's lemma to make higher order approximations or by local linearisation of the coef- 
ficients ll24l [301 ). but all suffer from common problems. In particular, minimising the approximation error (by 
increasing the mesh density) comes at the expense of increased computational cost, and further approximation or 
interpolation is needed to obtain the sample path at non-mesh points (which can be non-trivial). As a consequence, 
simulating sample paths by discretisation can be ill suited for subsequent application to a particular task. 

Recently, a new class of Exact Algorithms for simulating sample paths at finite collections of time points without 
approximation error have been developed for both diffiisions ll9l l6l l7l [T3l and jump diffusions lfT2l ITtI l20l . These 
algorithms are based on rejection sampling, noting that sample paths can be drawn from the (target) measure T'' 
by instead drawing sample paths from an equivalent proposal measure P'', and accepting or rejecting them with 
probability given by the Radon-Nikodym derivative of T'' with respect to P''. However, as with discretisation 
schemes, given a simulated sample path at a finite collection of time points subsequent simulation of the sample 
path at any other intermediary point may require approximation or interpolation and may not be exact. 

In this paper we introduce a novel mathematical framework for constructing exact algorithms which address this 
problem. In particular, instead of exactly simulating sample paths at finite collections of time points we instead 
focus on the extended notion of simulating skeletons which in addition characterise the entire sample path. 

Definition 1.1 (Skeleton). A skeleton (S) is a finite dimensional representation of a diffusion sample path (V ~ 
T' j, simulated without any approximation error by means of a proposal sample path drawn from an equivalent 
proposal measure fP '' j and accepted with probability proportional to -^^7, which is sufficient to restore the sample 
path at any finite collection of time points exactly with finite computation where V{cLi)\S ~ P''|<S. A skeleton 
typically comprises information regarding the sample path at a finite collection of time points and path space 
information which ensures the sample path is almost surely constrained to some compact interval. 

Methodology for simulating skeletons (the size and structure of which is dependent on exogenous randomness) is 
driven by both computational and mathematical considerations (i.e. we need to ensure the required computation 
is finite and the skeleton is exact). Central to both notions is that the path space of the proposal measure P'' can 
be broken into a discrete partition (a set of layers), and that the layer any particular sample path belongs to can be 
simulated. 

Definition 1.2 (Layer). A layer R(V), is is a function of a diffusion sample path V ~ P'' which determines the 
compact interval to which any particular sample path V is almost surely constrained. 

We show that, unlike existing exact algorithms, a valid exact algorithm can be constructed if it is possible to 
partition the proposal path space into layers, simulate unbiasedly which layer a proposal sample path belongs to 
and then, conditional on the layer, simulate a skeleton. Our exact algorithm framework for simulating skeletons is 
based on three principles for choosing a proposal measure and simulating a path space layer, 

PI - Principle 1 (Layer Construction). The path space of the proposal measure P'', can be discretely partitioned 
and the layer a proposal sample path belongs to (i.e. R{V) ~ "R) can be unbiasedly simulated. 

P2 - Principle 2 (Proposal Exactness). We can simulate from the probability law P'' o y^' o /« particular, 
conditional on Vq — v, Vj and R(V), we cam simulate any finite collection of intermediate points of the 
trajectory of the proposal diffusion exactly. 

Together PI and P2 ensure it is possible to simulate a skeleton. However, in addition we want to characterise the 
entire sample path and so we construct exact algorithms with the additional principle P3. 

P3 - Principle 3 (Path Restoration). Any finite collection of intermediate (inference) points, conditional on the 
skeleton, can be simulated exactly (i.e. we can sample Vt^{(jS), . . . , Vf„(<y) ~ P^|.Sj 

We also introduce a novel class of Adaptive Exact Algorithms for diffusions and jump diffusions, supported by 
new results for the simulation of Brownian path space probabilities (which are of separate interest) and layered 
Brownian motion (Brownian motion conditioned to remain in a layer). In addition, these results allow us to make 
methodological improvements to both the existing and new class of exact algorithm. 

We apply the results developed in this paper to significantly extend e-Strong Simulation methodology fSl (which 
allows the simulation of upper and lower bounding processes which almost surely constrain stochastic process 
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sample paths to any specified tolerance), from Brownian motion sample paths to a general class of jump dif- 
fusions, and introduce novel results to ensure the methodology in [8 | can be implemented exactly. Finally, we 
provide a number of examples along with simulation results. 

In summary, the main contributions of this paper are as follows: 

• A mathematical framework for constructing exact algorithms and simulating diffusion and jump diffusion 
sample path skeletons (see Sections [3] and |4]i, under weaker conditions (see Sections |2]and|3]l. 



An extension of existing exact algorithms to satisfy P3 (see Sections 3.1 4. 1 4.2 and|6]l, including a number 



of methodological improvements (see Sections [3| |3 . 1 1 |4T| |4. 3 [ |4i4| and 5.2 1 



A new class of adaptive exact algorithms for diffusions and jump diffusions (see Sections 3.2 4.3 and|7]l. 



• Methodology for simulating unbiasedly events of probability corresponding to various Brownian path space 
probabilities (see Section[5]l. 

• Methodology for the e-strong simulation of Brownian motion sample paths to ensure it can be initialised 
and is exact (see Sections [7]and[8]l. 

• An extension to e-strong simulation methodology from Brownian motion sample paths to a general class of 
jump diffusions (see SectionjH]). 

This paper is organised as follows: In Section [2] we detail conditions sufficient to establish results necessary 
for applying the methodology in this paper. In Sections [3] and |4] we outline our exact algorithm framework for 
diffusions and jump diffusions respectively, including the adaptive exact algorithms. We extend existing layered 
Brownian bridge constructions in Section |6] introducing novel constructions for the adaptive exact algorithms in 
Section [T] (both of which rely on novel Brownian path space simulation results which are summarised in Section 
|5]l. Finally, in Section[8]we apply our methodology to the e-strong simulation of (jump) diffusions. 

2 Preliminaries & Conditions 

The results underpinning this paper are for diffusions that are solutions to SDKs with unit volatility coefficient, 

dX, = a{X,) dt + dW„ = x € R, te[Q,T], (3) 
where W, is Brownian Motion and a: R ^ IR is the drift function and satisfies the following (mild) conditions, 
CI - Condition 1. a has a linear growth bound, i.e. 3K > Q such that \a{x)\ < ^'(l + |ji:|). 
C2 - Condition 2. a is continuously dijferentiable, i.e. a € C'. 

Applying the Mean Value Theorem, C2 ensures a is locally Lipschitz and so (|3]l admits a unique weak solution 
||27| . Allowing to denote the measure induced by pi and as Wiener measure (on the interval [0, T] with 
Xq - x) and A(m) := J^^ Q'(3')d3', we have, 

Rl - As a result of CI, the Radon-Nikodym derivative of Q ' with respect to exists and is given by Girsanov's 
formula (see for instance Il5l l27ll ). Furthermore, as a consequence of C2, we have A e and so we can 
apply Ito's Lemma to remove the stochastic integral (denoting (piX,) := a^(Xs)/2 + a'(Xs)/2), 



d(Q' 
dW-' 



(X) = exp|j^ a(X,)dX,-^ ^c?(X,) ds^ ^ exp i^MXr) - A(x) - 0(X,)dij 



(4) 



R2 - As a consequence of CI we have that A has a quadratic growth bound and so there exists some Tq < co such 
that Vr < To: 

c(y; X, To):^ exp i^A(y) - ^^^^| dy < oo. (5) 

R3 - By C2, a and a' are bounded on compact sets. In particular, suppose 3{,u e R such that V f e [0, T], 
X,(a)) e [£, v] 3Lx :^L (X(w)) € R, f/^ := f/ {X(oj)) e R such that V f e [0, T], (X,(w)) 6 [Lx, Uxl 
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C1-C2 are sufficient to ensure that R1-R3 hold (so we can implement the methodology in this paper), but can be 
relaxed if the Radon-Nikodym derivative of with respect to a proposal measure exists and can be bounded (for 
instance see |fT6ll ). 



Now, in order to tackle jump diffusions and heteroskedastic processes (as per ([T]i), further conditions are imposed, 

C3- Condition 3. yS € C'. 

C4 - Condition 4. cr e and strictly positive. 

C5 - Condition 5. A is locally bounded. 

Together with conditions CI and C2, C3-C5 are sufficient for the Lamperti transformation for jump diffusions 
Il24l [30l to be employed as follows (allowing methodology for the simulation of diffusions with unit volatility (j3]l 
to be extended to jump diffusions), 

R4 - Let rjiVi) -: X, be a transformed process, where ri(Vt) :- J^,, l/cr(u)du (where v* is an arbitrary element in 
the state space of V). Denoting by A^, :- 2,>i 1,,<, a Poisson jump counting measure (with respect to ff') 
and applying Ito's formula for jump diffusions to find dX, we have, 

dx, = [jj' dv; + 7j"{dv;f /2] + [jjiVr. + v(v,-)) - i](v,.)] dN, 

dt + dWr + (77 [r]-' (X,^ + v{t]-' (X,^))] - X,.) dNr . (6) 



This transformation is typically possible for univariate diffusions and for many multivariate diffusions [1]. A sig- 
nificant class of multivariate diffusions can be simulated by direct application of our methodology, and ongoing 
work is aimed at extending these methodologies more broadly to multivariate diffusions (see |34|). 

Denoting by and W \ the measures induced by (j6| and a driftless version of (|6| respectively, allowing 
4^1, ... , iJ/Nj. to denote the jump times in the interval [0, T], tfro := and i//Nt+i- ■- 4'Nt+\ '■- T, we have that the 
following holds: 

R5 - Under C1-C2 the Radon-Nikodym derivative of Q" with respect to W exists and is given by Girsanov's 
formula for jump diffusions 1,28. i24J . 



_dQ_'-* 
dW' 



-(X) = exp f^A(XT) - A(x) - cpiX,) d^ - g [A(X^,) - A(X^,_)l| . (7) 



3 Exact Algorithms 

In this section we outline how to simulate skeletons for diffusion sample paths which can be represented (under 
the conditions in Section [2] and following the transformation in (|6]l), as the solution to SDEs with unit volatility, 

dX, = a{X,) dt + dW„ = x € R, te[Q,T]. (8) 



As discussed in Section [T] we can view exact algorithms as a class of rejection samplers operating on diffusion 
path space. In this section we begin by introducing rejection sampling and outline an (idealised) rejection sampler 
for simulating entire diffusion sample paths. However, for computational reasons, this idealised rejection sampler 
can't be implemented so we instead, with the aid of new results and algorithmic step reordering, address this issue 
and construct a rejection sampler for simulating sample path skeletons which only requires finite computation (a 
retrospective rejection sampler). We present two algorithmic interpretations of this retrospective rejection sampler 
In Section 3.1 we present the Bounded and Unbounded Exact Algorithms, which are methodological extensions 
of existing exact algorithms ifTl [TtII . Finally, in Section 3.2 we introduce the novel Adaptive Unbounded Exact 
Algorithm. 
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Rejection sampling Il33l is a standard Monte Carlo method in which we can sample from some (inaccessible) 
distribution n, provided there exists an accessible distribution g with respect to which n is absolutely continuous 
with bounded Radon-Nikodym derivative. In particular, if we can find a bound M such that sup^ -g^ < M < oo, 
then drawing X ~ g and accepting (/ = 1) the draw with probability PgiX) ^ -g(X) € [0, 1] then {X\I = 1) ~ tt. 

An (idealised) rejection sampler to draw entire sample paths from (our target measure, induced by ([8])) can be 
analogously constructed using proposal sample paths simulated from Wiener measure W (which is the natural 
equivalent measure as (|8]l has unit volatility) which are accepted with probability Py^4.X) (which is proportional 
to the Radon-Nikodym derivative of (Q^ with respect to as given in (j4|, noting that A{x) is simply a constant). 

Idealised Rejection Sampler 



1 - Simulate X ~ W' 

2 - With probability PwK^^) = 6xp < A(X7-) - J^^ <p{Xs) ds \ /M accept, else reject and return to 1 



Unfortunately, this idealised rejection sampler can't be implemented as it isn't possible to draw entire sample paths 
from W (they're infinite dimensional random variables). Furthermore, (p{Xs) is typically unbounded and A only 
has a quadratic growth bound (see R2), so an appropriate bound (M < oo) to evaluate the acceptance probability 
can't typically be found. Finally, it is not possible to evaluate the integral expression in the acceptance probability. 

Clearly, as we can't simulate entire sample paths from Q^, the best we can hope for is to simulate a sample 
path skeleton. A retrospective rejection sampler (which only requires finite computation) can be constructed if it 
is possible to evaluate the acceptance probability of a proposal sample path without generating it in its entirety. If 
this were possible we could instead simulate a finite dimensional subset of the proposal sample path (a skeleton) 
then, given it is accepted, simulate any other aspect of the sample path required. The key to achieving this is 
noting that we don't need to evaluate the acceptance probability (which would require the entire sample path), but 
instead it is sufficient to find and simulate a binary random variable /, such that P(/ - l\S) - P^AX). 

We begin by finding a suitable bound (M < oo) for the acceptance probability. To remove the unbounded function 
A{Xt) from the acceptance probability we first simulate the end point from the law of Biased Brownian Motion 
(BBM) |9| and consider the resulting modification to the acceptance probability. 

Definition 3.1 (|9, Prop. 3]). Biased Brownian motion is the process Z, = (W,|Wo — x, Wt :—y~h) with measure 
TL^ , where jc, y e R, f e [0, T] and h is defined as (by R2 VT < < oo, h(y; x, T) is integrable), 

h(y;x,T):^ —^—exp{A(y)-^^^\, (9) 
c(y;x,T) [ IT j 

Theorem 3.1 (Biased Brownian Motion |9 , Prop. 3]). Q'* is equivalent to W"* with Radon-Nikodym derivative: 



dQ\ 



T 



^^-(X)(xexp^- ( m,)ds}. (10) 

Sample paths can be drawn from in two steps by first simulating the end point Xj :- y ~ h (although h doesn't 
have a tractable form a rejection sampler with Gaussian proposal can typically be constructed), then simulating the 
remainder of the sample path in the interval (0, T) from the law of a Brownian bridge, (X\Xq - x,Xt -y) ~ W^'^. 



As a direct consequence of ( 10 1, the acceptance probability of sample paths proposed from W is as follows. 



Pzv(X) = exp^- j| 4,{X,;)ds) IM. (11) 

Now, to determine whether a simulated proposal sample path is accepted we construct and simulate a binary 
random variable /, such that P(/ = 1|>S) = P%A.X). This can be achieved by means of an unbiased estimator of 
P'Z'iX) which is almost surely constrained to the interval [0, 1]. Recalling that 0(X[o,7-]) is bounded on compact 
sets (as a consequence of R3) note that if (after simulating the end point from BBM) we partition the path space 
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of 'W\Xj into disjoint layers and unbiasedly simulate which layer our proposal sample path belongs to (as in PI, 
denoting Rx ~ the simulated layer), then an upper and lower bound for 0(X[o,7-]) can be found conditional 
on this layer {Ux 6 R and Lx £ R respectively). Now, letting Ks(x) be the law of a- ~ Poi((f/x - Lx)T), U^^ the 
distribution of (^i, . . . ~ U[Q, TY, and considering a Taylor series expansion of the exponential function in 
( [TT| i we have. 



= E. 



Ux - (piX,) 
{Ux - Lx)T 



exp {-{Ux - Lx)T] [{Ux - Lx)T]' 

n 



ds 



IM 



= Ej? 



"^^^EK;;(;()EtU, 



(X {Ux-LxW'^''^ 
Ux - 0(^ft) 



U^ 



IM 



(12) 
(13) 



Now, as we only need to find a bound M < oo such that our unbiased estimator of P'Z'{X) lies in [0, 1], then we 
can choose M :- e^^"^ . The proposal sample path can be accepted or rejected with probability given by ( 13 1. 
Details on how to simulate layer information and simulate intermediary points conditional on layer information is 
presented in Sections [6]and[7] Note that in order to evaluate ( 13 1 we only need to simulate layer information and k 
{< oo) intermediary points, so we can now construct a retrospective rejection sampler (exact algorithm) requiring 



only finite computation. We present two distinct interpretations of this. In Section 3. 1 we present the Bounded and 



Unbounded Exact Algorithms and in Section 3.2 we present the Adaptive Unbounded Exact Algorithm. It should 



be noted that as a consequence of ( 13 i we don't require an explicit lower bound for (p{X[Qj]) (unlike existing exact 
algorithms ||9l|6l|7][l2l|20l), instead we note that under the conditions in Section [2] there almost surely exists a 
layer dependent lower bound which suffices. 

3.1 Bounded & Unbounded Exact Algorithms 

The Unbounded Exact Algorithm presented in Algorithm 3.1 is a retrospective rejection sampler directly based 



on simulating a finite dimensional proposal skeleton as suggested by ( 13 1, accepting the proposal with probability 



given by ( 13 1 and afterwards infilling the sample path at any other desired finite collection of time points. 



Algorithm 3.1 - Unbounded Exact Algorithm (UEA) 



1 - Simulate skeleton end point Xj : — y ~ h 

2 - Simulate layer information Rx ~ "R 

3 - Simulate skeleton points (x^^ , ■ ■ ■ , X^^^ \Rx 

3.1- Simulate k ~ Poi((f/x — l^x)T) and skeleton times , . . . , ~ U\Q, T] 
3.2 - Simulate sample path at skeleton times X^^ , . . . , X^^ ~ W^'' \Rx 

4 - With probability Wi=\ [(^-f ~ 4'{^^i)) I (Ux — ^x)] , accept entire path, else reject and return to 1 



* - Simulate inference times X,^,. . . , X,^^ ~ ^®^^/ W^^'^ j \Rx 



A number of alternate methods for simulating unbiasedly layer information (Step 2), layered Brownian bridges 
(Step 3), and the sample path at further inference times after acceptance (Step *), are given in Section|6] The skele- 
ton of the accepted sample path includes any information simulated for the purpose of evaluating the acceptance 
probability (any subsequent simulation must be consistent with the skeleton). As such, the skeleton is composed 
of the start and simulated end point {X^„ :- x and X^^^, :- y), skeletal points {X^, , . . . , X^J and layer Rx, 

SvEA {X{cj)) : = {(^,-, X^, J.^J^ ,RxY (14) 

The UEA extends the notion of an exact algorithm beyond simulating a sample path at a finite collections of time 
points without approximation, to the simulation of skeletons which characterise the entire sample path and can be 
used to restore the sample path at any finite collection of inference time points. Any valid combination of layer and 
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layered Brownian bridge can be inserted in Algorithm 3. 1 to reconstruct any existing exact algorithm for diffusions 
lIH |6l 111 [El (which satisfy P1-P2). However, unlike existing exact algorithms, the UEA in addition satisfies P3. 
Step * of Algorithm 3.1 can't be conducted in existing exact algorithms as the layer (Rx) imparts information 
across the entire interval. In Section [6] we show that Step * is possible (with some additional computation), by 
augmenting the skeleton with sub-interval layer information (denoting R^"''^^ as the layer for the sub-interval [a, b]). 



(15) 



The augmented skeleton allows the sample path to be decomposed into conditionally independent paths between 
each of the skeletal points and so the layer Rx no longer imparts information across the entire interval [0, T]. As 
such additional inference points required in Step * of Algorithm 3.1 can now be directly simulated. 



(16) 



In the particular case in which 0(X[o j]) is almost surely bounded there is no need to simulate layer information in 
Algorithm 3.1, the skeleton ( [T7| can be simulated from the law of a Brownian bridge and given the skeleton we 
can directly simulate inference points (so we satisfy P3 and no skeleton augmentation is required). This leads to 
the Exact Algorithm 1 (EAl) originally proposed in |6|, which we term the Bounded Exact Algorithm (BEA). 



>SBEA(X(w)):={(^i,Xft)^^'^3'}. 



(17) 



Diffusion sample paths can be simulated by successive simulation of sample paths of shorter length over the 
required interval by applying the Markov property, noting the Radon-Nikodym derivative in Q decomposes as 
follows (for any t e [0, T]), 



dQ^ 

dW^ 



(X) = 



exp\A{X,)-A(x) 



Jo 



^*(X,)ds 



exp M(Xr) - A(X, 



0(X,)di 



(18) 



A direct application of this decomposition is used in Section]?] As the computational cost of simulating a sample 
path doesn't scale linearly with interval length this decomposition is widely applied in practise. 



3.2 Adaptive Unbounded Exact Algorithm 

Within this section we outline a nowsl Adaptive Unbounded Exact Algorithm (AUEA). To motivate this we revisit 
( 12|13 1 and consider the Unbounded Exact Algorithm after the simulation of layer information. 



< exp {- {Ux - Lx) T] < exp |- j| [0(X,) - Lx] dij = Pz>|R(X)(^) < 1. 



(19) 



Reinterpreting the estimator in ( 13 1, we are exploiting the fact that P7i,'\r{x){^) is precisely the probability a Pois- 
son process of rate 1 on the graph Qa '■- {{a, b) : a e [0, T], be [Q,(p (^[oj]) - Lx]} contains no points. As this is 
a difficult space in which to simulate a Poisson process, we are instead simulating a Poisson process of rate 1 on 
the larger graph 0p := {(a, b) : a e [Q,T], b e [0, Ux - Lx]} (which is easier as Ux - Lx is a constant) and then 
conducting Poisson thinning (accepting the entire sample path if there are no Poisson points in 0^ £ &p)- 



As 0p can be much larger than 0a the resulting exact algorithm can be inefficient and computationally expensive. 
In this section we propose an adaptive scheme which exploits the simulation of intermediate skeletal points of the 
proposal sample path in Algorithm 3.1. In particular, note that each simulated skeletal point implicitly provides 
information regarding the layer the sample path is contained within in both the sub-interval before and after it. As 
such, by simulating each point separately we can use this information to construct a modified proposal c ^p, 
composed of a Poisson process with piecewise constant intensity, for the simulation of the remaining points. 

In Algorithm 3.1, the simulation of a proposal Poisson process of intensity Ax on the interval [Q,T] to deter- 
mine the skeletal time points (^i , . . . , ^k), can alternatively be simulated by exploiting the exponential waiting 
time property between successive events |23|. In particular, denoting Ti, . . . , as the time between each event 
^1, . . . then the timing of the events can be simulated by successive Exp(Ax) waiting times while 2; Tj < T. 
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The independence of the points of a Poisson process allows us to simulate them in any convenient order. In 
our case it is likely the sample path at points closer to the mid-point of the interval will contain more information 
about the layer the entire sample path is contained within. As such, there is an advantage in simulating these points 
first. If we begin at the interval mid-point (T/2), we can find the point closest to it by simulating Exp(2Ax) random 
variable t (we are simulating the first point at either side of the mid-point). As such, the simulated point (denoted 



^) will be with equal probability at either T jl - t ox T jl + t. Considering this in the context of ( 19 1, upon simu- 
lating ^ we have simply broken the acceptance probability into sub-intervals, the realisation of the sample path at 
the simulated time (X^) providing a binary unbiased estimate of the central sub-interval (where u ~ U[Q, 1]), 



Pzmx,(X) = exp |- j| ' [cpiXs) - Lx] ds - j^^^ [0(X,) - Lx] dsj E {l 



u < 



Ux - <^(Xf ) 



(20) 



If the central sub-interval is rejected the entire sample path can be discarded. However, if it is accepted then the 
acceptance of the entire sample path is conditional on the acceptance of both the left and right hand sub-intervals 



in ( 20 1, each of which has the same structural form as we originally had in ( 19 1. As such, for each we can simply 



iterate the above process until we have exhausted the entire interval [0, T]. 



As outlined above our approach is an algorithmic reinterpretation, but otherwise identical, to Algorithm 3.1. How- 
ever, we now have the flexibility to exploit the simulated skeletal point X^, to simulate new layer information for 
the remaining sub-intervals conditional on the exis ting layer Rx (as outlined in more detail in Sectionj7|. In par- 
ticular, if we evaluate the left hand sub-interval in (20 1, we can find new layer information (denoted Ry^^) which 
will contain tighter bound information regarding the sample path (fx < < X[()_^](w) < v^^'^ < vx) and so (as 
a consequence of R3) can be used to infer tighter bounds for (^(X[o,f]) (denoted f/^'^^ (< Ux) and i^''^^ (> Lx)), 



^'°'^\„,i (X) = expi- 



^[0(X,)-Lx]cl^Uexp|-(4°'^l-Lx)(|-r)}exp|- P [^(X,)-^''^] d* I . (21) 



The left hand hand exponential in (21 1 is a constant and it is trivial immediately reject the entire path with the 



complement of this probability. The right hand exponential of (21 1 has the same form as ( 19 1 and so the same 



approach as outlined above can be employed, but over the shorter interval [0, T jl — t] and with the lower rate 
(< Ax). As a consequence, the expected number of intermediary points required in order to evaluate the 



Am 



acceptance probability in ( 19 1 is lower than the UEA in Algorithm 3.1. 



This leads to the novel Adaptive Unbounded Exact Algorithm (AUEA) detailed in Algorithm 3.2. We outline 
how to simulate (unbiasedly) layer information (Step 2), intermediate skeletal points (Step 3.1.2) and new layer 
information (Step 3.1.4) in a variety of ways in Section [7] Our iterative scheme outputs a skeleton comprising 
skeletal points and layer information for the intervals between consecutive skeletal points. The AUEA with this 
skeleton has the distinct advantage that P1-P3 are satisfied directly. In particular, any intermediate inference 
points required after the skeleton has been simulated can be simulated directly (by application of Step 3.1.2 and 
Step 3.1.4 in Algorithm 3.2), without any augmentation of the skeleton (as in Algorithm 3.1). 



^auea(^(^^)) := ,«^'-''^''); 



(22) 



In Algorithm 3.2 we introduce simplifying notation, motivated by its recursive nature in which (as shown in (20 1) 



the acceptance probability is iteratively decomposed into intervals which have been estimated and are yet to be 
estimated. 11 denotes the set comprising information required to evaluate the acceptance probability for each of 
the intervals still to be estimated, 11 :- {YVi]^^^ Each 11, contains information regarding the time interval it applies 
to, the sample path at known points at either side of this interval and the associated simulated layer information, 
which we denote [sn(,> fn(;)], xn(o(:= ^"i''), 3'n(o(:= ^,+''') and respectively (where s^^^^ < sn(o < fn(/) < fn(,))- 
?n(0 ,,caH A\r-ar.tu, \r^far- >,rM,r.Hc f^^ A. fQj- (jjig speclfic Sample patfi over the interval [sncO' ^nco] 

= (*n(o + fn(o) /2, dn(i) := (fn(0 - *n(o) /2 and S Hi. 



As before, R^ can be used to directly infer bounds for ( 



(namely L^'^K 



f/J*'^ and A^*'^). We further denote mn(,) : 
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Algorithm 3.2 - Adaptive Unbounded Exact Algorithm (AUEA) 



1 - Simulate skeleton end point Xj : — y ~ h 

2 - Simulate initial layer information Rx ~ "R, setting H :- {E) :- {{[0, T],Xq,Xt,Rx}} and k - 

3 - While |n| 9t 0, 

3.1 - Simulate r ~ Exp(2A^). If t > d^ then set 11 := 11 \ H else 

3.1.1 - Set K — K + I and with probability 1/2 set £^'^ — — t else = + t 
3.1.2- Simulate Xe ~ W'^^iff' It?! 

3.1.3 - With probability {\ — |^J7^ — (p (^f^)] /^x) ^^j^ct path and return to 1 

3.1.4- Simulate new layer information axid f^^"'^^"^^ conditional on R^ 

3.1.5 - With probability {\ — exp |— l}^ ^. ^fc-'+^-'l _ _ reject path and return to 1 

3.1.6 - Set n := nu |[iH,'«H-T] (J |['MH+T,fH] ,%,^f+,/?[^'''^'^'"'} \ ^ 

4 - Define skeletal points . . . in temporal order from the set j^J , . . . , 
* - Simulate inference times X,, , . . . , X,„ ~ {^'^^^W^''-''^^' \Rf'''^'^^ 



4 Jump Exact Algorithms 



In this section we extend the methodology of Section |3] constructing exact algorithms for simulating skeletons of 
jump diffusion sample paths which can be represented as the solution to the following SDE, 



dX, = a (X,.) dt + dW, + djf''' Xo^xeM., te[0, T] 



(23) 



Denoting as the measure induced by (23 i, we can draw sample paths from \ by instead drawing sample 



paths from an equivalent proposal measure (a natural choice being that induced by the driftless SDE in (24 1), 
and accepting them with probability proportional to the Radon-Nikodym derivative of Q" with respect to W'^. 



dX, = dW, + dJ, '' = X € R; te [0, T] 



(24) 



The resulting Radon-Nikodym derivative (j7| differs from that for diffusions (|4|i with the inclusion of an additional 
term, so the methodology of Section[3]can't be directly applied. However, (j7|i can be re-expressed in a product form 
similar to ( 4| 1 8 1 (with iff i, ... , ij/M^ denoting the jump times in the interval [0, T], with i/tq ;= and 4'Nt+\ '■- T), 



dQ" 
dW" 



(X) 



Nt + I 

n 



expM(X^,_)-A(X^,_,)- 



<P(X,) ds 



(25) 



This product form of the Radon-Nikodym derivative is the basis for constructing Jump Exact Algorithms. Recall 
that decomposing the Radon-Nikodym derivative for diffusions justified the simulation of sample paths by succes- 
sive simulation of sample paths of shorter length over the required interval (see ([T8]l). Similarly, jump diffusion 
sample paths can be simulated by simulating diffusion sample paths of shorter length between consecutive jumps. 



In this section we present three novel Jump Exact Algorithms (JEAs). In contrast with existing algorithms 
lfT2l [TTl 120 1). we note that the Bounded, Unbounded and Adaptive Unbounded Exact Algorithms in Section [3] 
can all be incorporated (with an appropriate choice of layered Brownian bridge construction) within any of the 



JEAs we develop. In Section 4. 1 we present the Bounded Jump Exact Algorithm, which is a reinterpretation and 
methodological extension of 1 12 1, addressing the case where there exists an explicit bound for the intensity of the 
jump process. In Section [4T| we present the Unbounded Jump Exact Algorithm which is an extension to existing 
exact algorithms IIT7ir20l in which the jump intensity is only locally bounded. Finally, in Section 4.3 we introduce 



an entirely novel Ac/a/7f;Ve Unbounded Jump Exact Algorithm based on the adaptive approach of Section 3.2 
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4.1 Bounded Jump Exact Algorithm 



The case where the jump diffusion we want to simulate d23jl has an explicit jump intensity bound (sup^gjo j-j A(Xi,) < 
A < oo) is of specific interest as the proposal jump times can be simulated in advance. In particular, proposal jump 
times, ^Fi, . . . , ^F^A can be simulated according to a Poisson process with the homogeneous bounding intensity A 
over the interval [0, T]. A simple Poisson thinning argument Il23l can be used to accept proposal jump times with 
probability A(X^i.)/A. As noted in 1 12|, this approach allows the construction of a highly efficient algorithmic in- 
terpretation of the decomposition in (25 i. The interval can be broken into segments corresponding precisely to the 
intervals between proposal jump times, then iteratively between successive times, an exact algorithm (as outlined 
in Section[3]l can be used to simulate a diffusion sample path skeleton. The terminal point of each skeleton can be 
used to determine whether the proposal jump time is accepted (and if so a jump simulated). 



The Bounded Jump Exact Algorithm (BJEA) we outline in Algorithm 4.1 is a modification of that originally 
proposed |12 | (where we define ^Fq := and "Paja+i := T). In particular, we simulate the proposal jump times 
iteratively (exploiting the exponential waiting time property of Poisson processes f23l as in Section [3^, noting 
that the best proposal distribution may be different for each sub-interval. Furthermore, we note that any of the 
exact algorithms we introduce in Section |3]can be incorporated in the BJEA (and so the BJEA will satisfy at least 
P1-P2). In particular, the BJEA skeleton is a concatenation of exact algorithm skeletons for the intervals between 
each proposal jump time, so to satisfy P3 and simulate the sample path at further inference times (Step *), we 
either augment the skeleton if the exact algorithm chosen is the UEA (as discussed in Sectio ns 1 3. 1 1 and [6]), or, if 
the exact algorithm chosen is the AUEA then simulate them directly (as discussed in Sections |3.2| and[7]l, 

'Sbjea(W):= [J (26) 



Algorithm 4.1 - Bounded Jump Exact Algorithm (BJEA) 

1 - Set j = 0. While '¥j < T, 

1.1 - Simulate t ~ Exp(A). Set j - j + 1 and ^'j - *Py_i + t 

1.2 - Apply exact algorithm to the interval [^I*;-!, jAT^j, outputting skeleton »S^^ 

1.3 - // '¥j > T then set Xj = Xt- else 

1.3.1 - With probability AiXfi.)/ A set X^. := Xvp^_ + fy^X^fi^^^ else set Xv^. := Xij<,_ 

*- Simulate inference times X,,, ... ,Xt ~ ®'^^^ |®'^'t'w,*^''"'J iTJ^fi , ,| 



4.2 Unbounded Jump Exact Algorithm 

Considering the construction of a JEA under the weaker condition C5 (in which we assume only that the jump 
intensity in ( [23]l is locally bounded), it is not possible to first simulate the jump times as in Section pTT] However, 
(as in Sectionpland as noted in lfT7l|20l ). it is possible to simulate a layer /?x ~ ^> and then infer a jump intensity 
bound (A < Ax < oo) conditional on this layer. As such we can construct a JEA in this case by simply incorporat- 
ing the jump intensity bound simulation within the layer framework of the UEA and AUEA. 

The Unbounded Jump Exact Algorithm (UJEA) in Algorithm 4.2 is a JEA construction based on the UEA and 
extended from ||20] . The UJEA is necessarily more complicated than the BJEA as simulating a layer in the UEA 
requires first simulating an end point. As a result, although (conditional on a layer) there exists a Poisson process 
in which to simulate proposal jump times, we can no longer segment the interval into sub-intervals according to 
the length of time until the next jump. Instead we need to simulate a diffusion sample path skeleton over the entire 
interval (along with all proposal jump times) and then determine the time of the first accepted jump and simulate it 
(if any). If a jump is accepted another diffusion sample path has to be proposed from the last accepted jump time 
to the end of the interval. This process is then iterated until no further jumps are accepted. The resulting UJEA 
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satisfies P1-P2, however, as a consequence of the layer construction, the jump diffusion skeleton is composed of 
the entirety of each proposed diffusion sample path skeleton (we can't apply the strong Markov property to discard 
the sample path skeleton after an accepted jump). In particular, as the simulated layer covers the entire interval (the 
sub-interval before and after an accepted jump) it imparts information across the interval, so the skeletal points in 
the interval after an accepted jump contain information about the sample path before the accepted jump. 



The UJEA doesn't satisfy P3 unless the skeleton is augmented (as with the UEA outlined in Sections 3.1 and[6]), 
however, this is computationally expensive and is not recommended in practice (hence we have omitted Step * 
from Algorithm 4.2). Alternatively we could use the AUEA within the UJEA to directly satisfy P3, however it is 



more efficient in this case to implement the Adaptive Unbounded Jump Exact Algorithm in Section 4.3 



Algorithm 4.2 - Unbounded Jump Exact Algorithm (UJEA) 



1 - Set j — and — 

1.1 - Simulate skeleton end point Xj '.- y ~ h(y; X^ ., T — ipj) 
1.2- Simulate layer information Ri,, , ^, ~ % and infer Ai, , ^, 

1.3 - Simulate proposal jump times Nj''' ~ Poi (a^j^ j^(T - tl/j)^ and . . . , ^^^a.j ~ U[iffj, T] 

1 .4 - Simulate skeleton points and diffusion at proposal jump times ( X.j , . . . , X.i , X^^,j , . . . , X^j 1 iRi,, . 7,, 

1 .4. 1 - Simulate Kj ~ Poi ([t^^^,^ , j-j - L-l^^^ j ■ [T — tp j)^ and skeleton times , . . . , ~ U[4>j, T] 
1.4.2- Simulate sample path at times X.J, ... ,X.],X^^j, ... ,X^j ~ 'W')'\\Rv, , -n 

1.5 - With probability ^1 - n^^=i {^^xitp t] ~ ('^f^)) ^ (^^[(^ t] ~ ^x[<ij t\I ) > ^^j^^^ '^^d return to 1.1 
1.6- For i in 1 to Nj'^ 

1.6.1 - With probability /l(X^./)/A^j^^_^j set X^, = X^i, X^i := X^i_ + fy{x^i^, ipj+i ^P/, j ^ j + 1, 



and return to 1.1 



4.3 Adaptive Unbounded Jump Exact Algorithm 

The novel Adaptive Unbounded Jump Exact Algorithm (AUJEA) which we present in Algorithm 4.3 is based on 
the AUEA and a reinterpretation of the UJEA. Considering the UJEA, note that if we simulate diffusion sample 
path skeletons using the AUEA then, as the AUEA satisfies P3 directly, we can simulate proposal jump times 
afterwards. As such, we only need to simulate the next proposal jump time (as opposed to all of the jump times). 



which (as argued in Section 3.2 1, provides further information about the sample path. In particular, the proposal 
jump time necessarily lies between two existing skeletal times, ^_ < < so the layer information for that 
interval, R^' can be updated with layer information for each sub-interval R^ ' and R^ (the mechanism 
is detailed in Section |6]l. Furthermore, upon accepting a proposal jump time the sample path skeleton in the 
sub-interval after 'I' contains no information regarding the skeleton preceding *!* (so it can be discarded). As such, 
the AUJEA satisfies P1-P3 and the skeleton is composed of only the accepted segments of each AUEA skeleton, 

^AUJEA (X{oj)) y S^l^^^l'^ {X{u)) . (28) 
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Algorithm 4.3 - Adaptive Unbounded Jump Exact Algorithm (AUJEA) 



1 - Set j — and — 

2 - Apply AUEA to interval ^ij/j, rj, outputting skeleton S^^^'^j^ 

3 - Setk^O and = tf/j. While T^' < T, 

3.1 - Infer , 

3.2 - Simulate t ~ Exnl A'' 1. Set k - k + 1 and ^{ -^'^ , + t 

^\ X['Vl,T]j * k-i 

3.3 -If'¥l<T 

3.3.1 - Simulate X^^j ~ 



3.3.2 - Simulate ' andTJ^^' ana sef 5^^^^ := >S^^g^ U <X>j,/,/?^ '^x M 

3.3.3- With probability A{X^i) I . set X^j - X^j, X^i :- X^j_ + fylx^j], tj/j+x '■- ^i, retain 

k X[1'^_j,7'] k k k k \ kj 

^^AUEA^^ ' discard set j — j + \ and return to 2 

4 - Define skeletal points X\, ■ ■ ■ ,Xm in temporal order from the time points in .Saujea •= Ui^=i' "^auea^ 
* - Simulate inference times X,„ . . . , X,„ ~ t.'v^i)'^' Vx"'*'' 



4.4 An Extension to the Unbounded & Adaptive Unbounded Jump Exact Algorithms 

As we don't know the timing of the next jump in the UJEA and AUJEA, in contrast with the BJEA, we simulate 
diffusion sample paths which are longer than necessary (so computationally more expensive), then (wastefully) 
partially discard them. To avoid this problem we could break the interval into segments and iteratively simulate 



diffusion sample paths of shorter length over the interval (as in (18 i), thereby minimising the length of discarded 
segments beyond an accepted jump. However, the computational cost of simulating a sample path does not scale 
linearly with the interval it has to be simulated over, so the optimal length to decompose the interval is unknown. 

It is possible to extend the UJEA and AUJEA based on this decomposition and Poisson superposition li23l . In 
particular, if it is possible to find a lower bound for the jump intensity Ai € (0, A), then we can consider the target 
jump process as being the superposition of two jump processes (one of homogeneous intensity Ai and the other 
with inhomogeneous intensity A - A\,). As such we can simulate the timing of an accepted jump in the jump 
diffusion sample path under the homogeneous jump intensity Ai by means of a r ~ Exp(/}4,) random variable. If 
T e [0, r] then there is no need to simulate proposal diffusion skeletons over the entire interval [0, T], instead we 
can simulate them over [0,t]. Furthermore, we can modify the bounding jump intensity in the UJEA and AUJEA 
for generating the proposal jump times in the proposal diffusion sample path skeletons from Ax to Ax - Ai. 

5 Brownian Path Space Simulation 

In this section we present key results which we use to construct layered Brownian bridges in Sections |6] and 
[T] In Section 5.1 we outline a number of established results pertaining to the simulation of a variety of aspects 



of Brownian bridge sample paths. In Section 5.2 we consider known results (along with some methodological 



improvements) for simulating events corresponding to the probability that Brownian and Bessel bridge sample 



paths are contained within particular intervals. Finally, in Section 5.3 we present novel work in which we consider 



simulating probabilities corresponding to a more complex Brownian path space partitioning. Central to Sections 



5.2 and 5.3 are Theorem 5.1 and Corollaries 5.1 and 5.2 which together form the basis for simulating events of 



unknown probabilities which can represented as alternating Cauchy sequences. 
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Theorem 5.1 (Series Sampling IfMl Chapter 4.5]). An event of (unknown) probability p e [0, 1], where there 
exists monotonically increasing and decreasing sequences, (5^ : k € K>o) and (Sj, : k e K>o) respectively, such 
that lim^^oo i p and limi:^oo '\ p, can be simulated unbiasedly. In particular, a binary random variable 
P :— l(u < p) can be constructed and evaluated (where u ~ U[Q, noting that as there almost surely exists a 
finite K :- 'mf{k : u i (S^,S^)} we have l(u < p) - 1(m < S^). 

Corollary 5.1 (Linear Transformation). Probabilities which are linear transformations or ratios of a collection of 
probabilities, each of which have upper and lower convergent sequences can be simulated by extension of Theorem 
In particular, suppose f : R™ — > IR+ e C' such that | d// dui(u)\ > V 1 < i < m and u € R™ and that the 



5.1 



probability p :— f{p\, . . . ,/?,„) then defining the sequences (T,' : k e Wi>q) and fr,' : k € K>oj as follows. 



ju-^lS'f //d//d«, >o r.+ = i^^ if dfidui>Q 

« 15;'^ //d//dM, <0 ' * \ ifdflAui<0 ■ ^ ' 

we have that 5^ := f(T^'^ , . . . , T"™' ) is monotonically increasing and converges to p from below and 5^ :— 
f(T^'^, . . . , r™'^) is monotonically decreasing and converges to p from above. 

Corollary 5.2 (Retrospective Inversion Sampling Q Prop. 1]). If p can be represented as the limit of an alternat- 
ing Cauchy sequence (St ■ k € S>oj, then splitting the sequence into subsequences composed of the odd and even 
terms respectively, then each subsequence will converge to p, one of which will be monotonically decreasing and 



the other monotonically increasing, so events of probability p can be simulated by extension of Theorem 5.1 



Algorithm 5.1 outlines how to simulate events of probability p (applying Corollary 5.2 1, where p can be repre 



sented by an alternating Cauchy sequence (the even terms converging from below and the odd terms from above), 
= 5o < ^2 < 54 < 56 < ... < /5 < ... < 55 < 53 < 5i < 1. (30) 

Algorithm 5.1 - Retrospective Inversion Sampling 

1 - Simulate u ~ t/[0, 1] and set k — \ 
2- While ue{S2k,S2k+\),k^k+\ 
3 - If u < Sok then accept else reject 



5.1 Simulating Brownian Bridges and Related Processes 

The density of a Brownian bridge sample path W^f, at an intermediate time q e {s,t) is Gaussian with mean 
w :- X + {q - s){y - x)/{t - s) and variance cr^^, := [t - q)(q - s)/{t - s) (so can be simulated directly). The joint 
distribution of the minimum value reached by H^,'/, and the time at which it is attained (r, m), is given by Il22l . 

Tot- ^ A ^ A \»7 w ^ {w-x){w-y) j (w-xf (w-yf\ 

P(m e dw, T e dq\Ws = x,W, - y) oc exp < } aw dq. (31) 

V(f - q)Hq - s)^ I '2(q-s) 2{t-q)] 



Analogously the maximum (t, Wj :- m), can be considered by reflection. As shown in ||6l we can draw from (31 



as outlined in Algorithm 5.2 which extends ||6l noting that it is possible to condition the minimum to occur within 
a particular interval (i.e. simulate (r, m)\ (m e [ai, 02]) where ai < a2 < x A y). 

Conditional on a Brownian bridge sample path minimum (or maximum), the law of the remainder of the tra- 
jectory is that of a Bessel bridge, which can be simulated by means of a 3-dimensional Brownian bridge of unit 
length conditioned to start and end at zero as outlined in li3J and Algorithm 5.3 (maximum by reflection). 
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Algorithm 5.2 - Brownian Bridge Simulation at its Minimum Point (constrained to the interval [fli,fl2] where 
fl] < a2 < X A y and conditional onWs — x and W, - y (denoting IGau(jj, A) as the inverse Gaussian distribution 
with mean ju and shape parameter A) 

1- Simulate Ui ~ U [M{a2), M(ai)] where M (a) :- expi^-2^xy + ((xAy)-ay^ - ((xAy)-a) ■ (x+y)^/{t-s)^ 
and U2 ~ U[0, 1] 

2 - Set m := X - ^ ^Jiy - x)^ - 2(t - s) log(Mi) - (y - x)j /2 

X - m , I y - m (y - \ . 1 / x - m (x - 

3 - // M2 < r then V ~ IGau -, else IGau -, 

X + y — 2m \x — m t — s ) V \y — m t — 



ihy 

s 



. „ t - s 

4- SetT-.^ 



1 + V 



Algorithm 5.3 - (Minimum) Bessel Bridge Simulation (at time q G (s, /) given — x,Wi — y and Wj — in). 



1 - If q < T then r - s else r - t. Simulate b\ ,b2,bj, '~ N (o, ^ — — — ^ 



2- SetW,:^m+ ^/\7^\ ■ J^J^^L^^t^ + b]\ bj + bl 



5.2 Simulating Elementary Brownian Path Space Probabilities 

In this section we briefly outline results pertaining to the probability that a Brownian bridge sample path is con- 



tained within a particular layer l|2][3T| (Theorem 5.2 1 and how to simulate events of this probability |7| (Corollary 



5.3 1. Similarly in Theorems 5.3 and 5.4 we outline a result (first shown in Q) in which we show that the proba- 



bility that a Bessel bridge sample path is contained within a particular layer can be represented as an infinite series. 



Of particular importance for what follows is Corollary 5.5 in which we establish that it is possible (without 



any further assumption regarding the layer) to simulate events with a probability corresponding to the probability 



that a Bessel bridge sample path is contained within a particular layer by application Corollary 5.2 



Theorem 5.2 ( 1311 Theorem 3]). The probability that a Brownian bridge sample path W ~ Wj',\ remains in the 
interval [{, v] (i.e. V m e [s, t] Wu e [£, v]) can be represented as an infinite series, 

oo 

/;J(x, y):^¥(We[e,v])^l-Y, (J; y) - if^^if x, y)] , (32) 

where g^'^ij; x,y) := x,y) + g'^j'^^i j; -x, -y\ ip^;"(j, x,y) -.^ ^P^'ii j; x,y) + ^;J"''(;; -x, -y) and, 

f[f (7; x,y) exp [\v - e\j + {£ A v) - x) ■ [\v -{\j + {£Av)- y)| , (33) 

x,y) - exp (jy - {fj + \v-e\{x- y)"^ . (34) 

Corollary 5.3 (I?] Prop. 2]). yl'1{x,y) is an alternating Cauchy sequence, so events of probability y^' "(x,}') can 
be simulated by retrospective inversion sampling ( Corollary \5.2\ and Algorithm 5.1) using the following sequence, 

k 

Si, 1 - 2 0; ^^y) - 'P^:;(j^ x,y)} , S^^^ - S^ - s-f-" (^ + hx,y). (35) 

j=l 
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As shown in Q, Theorem 5.2 and Corollary 5.3 can be extended to consider simulating events with a probability 
corresponding to the probability a Bessel bridge sample path is contained within a particular layer As indicated 
in Definition |5.1 1 we have to consider two possible cases where either of the end points are attain the sample path 
minimum (or maximum) or not. 

Definition 5.1. We allow 6ff(x,y) := l{m < (xAy)} ■ 6ff(l;x,y) + l{m = (xAy)} ■ d'^fi2;x,y) to denote the 
probability that a Bessel bridge sample path W ~ W'J'^Iot, (with minimum in) remains in the interval [m, u] 



f V M e [i, f] W„ e [m, v\) where 6"^f{\; x,y') and 6'^f(2; x,y) are defined in Theorems 5.3 and 5.4 respectively. 



Remark 5.1. We can similarly consider the probability that a Bessel bridge sample path W ~ Wj',|m (with 
maximum m) remains in the interval [f, m] f V m e [s, t] Wu € [£, m]) by a simple reflection argument of the above. 

We first consider the case where the neither end point attains the Bessel bridge minimum (or maximum). 

Theorem 5.3 (|7., Prop. 3]). The probability that a Bessel bridge sample path W ~ W^','|m, (with minimum 
m < (xf\y)) remains in the interval [m, v] (W u € [s, t] Wu e [m, v]) can be represented as an infinite series, 

m,v/ \ 

^"(1; ^,3;) := Viw e [m, v]\w > m, (xAy) > m) = , rTTTT- (36) 

^ ' ' I - exp { - 2{x - m)(y -m)/t\ 

Corollary 5.4 (|7, Prop. 3]). Events of probability 6'^f(l; x,y) can be simulated by application of retrospective 
inversion sampling (as per Corollaries 5.1 5.2 and Algorithm 5.1) using the following sequence, 

Sf^ := ; — -. (37) 

1 - exp { - 2(x - m)(3; - m)/f) 

We now consider the case where the either of end points attains the Bessel bridge minimum (or maximum). 

Theorem 5.4 (|j7] Prop. 3]). The probability that a Bessel bridge sample path W ~ W'^'j^jm, (with minimum 
m — X < y) remains in the interval [m, y] (V m e [s, t] Wu e [m, v]) can be represented as an infinite series, 

6ff(2; x,y) P (w^ e [m, v]\W > m) = 1 - ^— ^ |] [r:f(fy) -/::r(i^y)] > (38) 

where we denote, 

Cr(i; y) -m\j-(y- m)) exp j ^_ ^ {jv-m\j-(y- m)j \ , (39) 

/II \ 1 - ml j . , I 

x^sfiL y) - y\y - '«| ; + (y - m)) exp j — 737^ (1'^ - "^1 ; + (}'- '«)) [ ■ (40) 

Remark 5.2 (fT, Prop. 3]). We can similarly consider the probability a Bessel bridge sample path W ~ W^f |m, 
(with minimum m — y < x) remains in the interval [m, v] by a simple refiection argument of Theorem\5 .4\ 



Corollary 5.5. After the inclusion of the first k :— sjit — s) + \v - mp/ (2\v — m\) terms, 5'^f{2; x,y) is an alter- 
nating Cauchy sequence, so events of probability 6"^'" (2; x,y) can be simulated by retrospective inversion sampling 



(as per Corollary 5.2 and Algorithm 5.1) using the following sequence (where A; e N such that k > k), 

- 1 - i [^f^j^y) -4rO;3')} ' ^2.1, - - + l;y). (41) 

Proof. As (y - m) G (0, {v - m)] then V j we have t//j,Xj ^ 0. As such it is sufficient to show that Vj>k that 
iffj > Xj ^ 'Pj+i ^ Xj+i ^ • ■ ■ which can be shown by first showing that Vj i/Zj > Xj and then '^jXj ^ 'Ay+i- 
Considering il^jlxj if j ^ ^ then this is minimised when y - m and if/jlxj > 1- Similarly considering Xj/^j+i if 
j >k then this is minimised when y - v where xj/ipj+i > 1 □ 
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Figure 



5.2 



y event -W e\l,v\ fh e \t, (xAy)], in e \{x\y), v\ 



5.3 Simulating Brownian Path Space Probabilities 

In this section we establish that the probability a Brownian bridge sample path, conditional on a number of in- 
termediate points, has a minimum in a lower layer and a maximum in an upper layer (or in each sub-interval a 
minimum in a lower layer and a maximum in an upper layer) can be represented as an infinite series (Theorems 



5.5 and 5.6 respectively), and events of this probability can be simulated (Corollaries 5.6 and 5.7 respectively). 



In this section we introduce the following simplifying notation, m,, := vc&\Wq,q e \s,t\\, ihs^ :- sup{W^;^ e 
[s,t]},Q {qi,...,q„}, ^ =wi,...,W,„ = w„}, £ e . . . ,m,„,, e [^^^^ ,]}, 

^ := [m.^qi e [v{g, , vl^^], ihq^^j e [v^^^,, v^,]], qo := s and q„+i := t. 

Theorem 5.5. The probability a Brownian bridge sample path W ~ W^f ["W has a minimum rhs^ e [^i, ^T] '^nd 
a maximum rhsj e [ui, ft] can be represented as an infinite series. 



(n) a/T,fi,fT 



{Q,^V) := P(m,,, e {a,t'[],m,j e [vi,vt]\-W) 



(42) 



=0 J L /=0 J L 1=0 J L ,=0 

Proof. Follows by sample path inclusion-exclusion and the Markov property for diffusions, 
'^"^P%];''^'"\Q,'\V) = P (iV e [a, v1]\w) - P (W e [(t uT]|"W) - P (w e [a, vUl^V) + ¥(We [ft vi]\w) 

n If" 



1=0 J L /=o 



□ 



Corollary 5.6. Events of probability p-"' can be simulated by retrospective inversion sampling ( as per Corallaries 



5.1 5.2 and Algorithm 5.1), noting that p^"^ function ofy probabilities, using the following sequence, 

n If" IT" 

j~j _jr(**+i.T.T) _ j~j^y(?„<?,+i.i.i) ^ j~j ^y(<?„*+i,T,J 



:P(") ._ 



1=0 



y(*,?,+i,i,T) 
k 



/=0 



1=0 



- i=0 



(43) 



Definition 5.2. We define p{s, q, t, x, w, y, £1, {t vi, l-T) := *'*p',w (<3' which is equal to p in 

Theorem 5.6. The probability a Brownian bridge sample path W ~ W,';' \'W, has in the sub-intervals induced 
between successive points in "W, a minimum and maximum in particular layers, £. and U respectively, can be 
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represented as an infinite series. 



/=() 

Proof. Follows by the Markov property for diffusions and sample path inclusion-exclusion, 

n 

!=() 

i=() 

Corolla ry 5.7. Events of probability yS*"' can be simulated by retrospective inversion sampling (as per Corollaries 
and Algorithm 5.1), noting that /3^"^ is a function ofy probabilities, using the following sequence. 



5.1 



5.2 



cP(n) ._ pr rny(9„?,+i,a,i'T) _ cr(9..*+l.^T,fT) _ gy(qi,qM,tLvi) _|_ ^yte,?,+ i/T,fi)] 
k ' Y k k+l k+\ k J" ^ ^ 



Deflnitloii 5.3. We define I3(s, t, x, y, ti, vi, u^) := *°^y8i;;",,((3, OV), which is equal to p in Ml/. 

6 Layered Brownian Bridge Constructions 

In this section we outline how to construct and simulate finite dimensional skeletons of layered Brownian bridges 
for use within the UEA (Algorithm 3.1), which is in turn used within the BJEA (Algorithm 4.1) and UJEA (Al- 
gorithm 4.2). In particular, we address how to simulate unbiasedly layer information (Algorithm 3. 1 Step 2), how 
to simulate intermediate skeletal points (Algorithm 3.1 Step 3) and how to simulate further inference times after 
acceptance of the proposed sample path (Algorithm 3.1 Step *). 

We present two alternate layered Brownian bridge constructions based on extensions to existing exact algorithms. 



In Section 6.1 we present the Bessel Approach, which is a reinterpretation of part of the Exact Algorithm 3 (EA3) 
proposed in |7|, in which we incorporate the methodological improvements outlined in Sections |3] and |5] and in- 
troduce a novel approach for conducting Algorithm 3.1 Step * (which could not previously be achieved). As a 
consequence, the resulting (complete) UEA, with the inclusion of the Bessel approach, satisfies P1-P3 (as op- 



posed to only P1-P2 in EA3 |7|). Finally, in Section 6.2 we outline a Localised Approach for constructing a 



layered Brownian bridge (based on the Localised Exact Algorithm (LEA) lfT3l \Tji), showing that the resulting 
UEA only satisfies P1-P2 and discussing the difficulties in conducting Algorithm 3.1 Step * and satisfying P3. 



In both the Bessel and Localised approaches it isn't possible to directly simulate intermediate points conditional 
on a simulated layer (as required in Algorithm 3.1 Step 2). Instead, proposal sample path skeletons are simulated 
by intermediate Monte Carlo techniques, including Rejection Sampling (see Section[3]l and demarginaUsation Il3?l . 

Demarginalisation ll33l is a technique whereby artificial extension of a density (with the incorporation of aux- 
iliary variables) simplifies sampling from it. To illustrate this consider the case where we want to draw a sample 
g{X), but this is not directly possible. However, suppose that with the introduction of an auxihary variable Y, 
sampling from g{X\Y) is possible and g{X, Y) admits g{X) as a marginal, 

g(X)^ r g{X\Y)g{Y)&Y (46) 

We can sample from giX) by first sampling Y from g{Y) and then sampling from g{X\Y). This algorithm can be 
viewed as a black box to generate samples from g{X) - Y can be simply marginalised out (i.e. 'thrown' away). 
Considering demarginalisation in the context of the exact algorithms, we can simulate any (auxiliary) aspect of 
the proposal diffusion sample path in addition to the skeleton to aid sampling. Provided the auxiliary information 
does not influence the acceptance probability then it is not part of the skeleton and doesn't need to be retained. 
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6.1 Bessel Approach (EA3) 



The central idea in the Bessel Approach is that proposal Brownian bridge sample paths can be simulated jointly 
with information regarding the interval in which they are contained (Algorithm 3.1 Step 2), by means of a par- 
titioning of Brownian bridge path space with an (arbitrary) increasing sequence, {fl,),>o, flo - 0, which radiates 
outwards from the interval [{x/\y), (xVy)] demarcating 'layers'. For instance, the i* layer is defined as follows. 



J, — {{x/\y) - Gi, (xVy) + fl,] . 



(47) 



The (smallest) layer, J = i, in which a particular Brownian bridge sample path is contained can be simulated by 
inversion and retrospective inversion sampling as in Algorithm 6.1. The CDF of I can be written as follows (with 
reference to Theorem |5. 2 [ and as shown in ||7|), 

(48) 



P(J < = ^{W^f e [(xAy) - fl„ (xWy) + a,]] ) = y(s, t,x,y, (xAy) - a,, (xWy) + a,). 



Algorithm 6.1 - Simulation of a Brownian Bridge Layer 



1 - Simulate u ~ U[0, 1] and set l — 1, A: = 

2 - While u e (Sl^j^^^ (s, t, x,y, (xAy) — ai, (xVy) + a,) , ^j^^ (s, f, x,y, (xAy) — at, (xVy) + a,)^, k — k + I 



3 - If u > SZ, set L — L + I and return to step 2 else set I — i and end 



Now, in order to simulate intermediate points (Algorithm 3.1 Step 3), we need to condition on the layer simulated 
in Algorithm 6. 1. In particular, denoting D, as the set of sample paths which are contained in the t* layer we have. 



where. 



(xVy), (xVy) + a, 



{xAy) - a„(xAy) - flj_i j| Q |m,_, 



(xAy) — fl,, (xAj 



(49) 

(50) 
(51) 



Directly simulating intermediate points from the law of D,, (denoted D,) is not possible. Instead (as proposed 
in 121) we can propose sample paths from the mixture measure IB, :- M,/2 + M,/2 (M, and M, being the law 
induced by M, and M, respectively) and accept them with probability given by the Radon-Nikodym derivative of 
ro, with respect to IB,, where. 



M, 



hs^t £ [(-^ Ay) — fl,, (xAy) - M, = \^s,t £ [(-"cVy) + (xVy) + a,j 



It was shown in fT\ that D, is absolutely continuous with respect to IB, with Radon-Nikodym derivative, 

dD, l(XeA) 

(^) ^ ^ 

dIB, 1 -H 1(X e M, n M,) 



(52) 



(53) 



Sample paths can be drawn from D, by proposing them from either M, or M, with probability 1 /2 and then 
accepting them with probability given by (53 i. For instance, conditional on simulating from M, we accept the 
sample path with probability 1 if the sample path maximum is contained within the (i- 1)* layer or with probability 
1 /2 if it is contained between the (i - 1)* and i* layer (and otherwise). To aid simulation and ensure we simulate 
from M,, we first simulate the sample path minimum Xj - fixs., as per Algorithm 5.2, and subsequently simulate 
any required skeletal points ^\,. . . from a Bessel bridge as per Algorithm 5.3. As we can only simulate our 



sample path at a finite collection of points we can't directly evaluate (53 i. However, we can find probabilities for 
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these events and so simulate whether or not they occur by application of Corollaries 5.3 and 5.4 and Lemmata 5.4 
and 5.5 (denoting X\^ - ■ ■ temporal order from the set {^i , . . . , ^i^, s, t, f)), 



K+2 



!=1 

K+2 

¥^{X e M, n M,) = F^^iX e A) - nC'i^i'"" (^^.'^^...) ' ^^^^ 



1=1 



As both ( 54 1 and ( 55 1 are probabilities which can be represented as a function of 6 probabilities, then events of 
this probability can be simulated by retrospective inversion sampling (as per Corollaries |5. 1 [ [S!2] and Algorithm 
5.1). The synthesis of the above approach for simulating Brownian bridge conditional on the layer simulated in 
Algorithm 6.1 (i.e. conducting Algorithm 3.1 Step 3) leads to Algorithm 6.2. 

Algorithm 6.2 - Layered Brownian Bridge Simulation (Bessel Approach) 

1 - Simulate U\,U2 ~ t/[0, 1], set j — k — 

2 - Simulate Auxiliary Information (conditional on I — l) 

2.1 - If ui < 1/2 simulate minimum Xj — fhs^, (— £i — £2) ond set Vi — (xWy) + and V2 — (xWy) + a, 

2.2 - If ui > 1/2 simulate maximum Xj — rhs^t i- u\ — V2) and set {\ — (xAy) — and I2 — (JcAy) — 

3 - Simulate intermediate times , , . . . , X^^ from a Bessel Bridge conditional on Xj 

4 - While U2 G ( n-:? s\.^i{iuv,), ntt 5^ .(a, f^i)), ; = ; + 1 

4.1 - If U2 < Y['il\ 'S2j+i{^i'^i) ^^^^ accept sample path 

4.2 -Ifu2> UiliStj{euVi) while U2 e ( 5^,^i(^2, ^2), Wi S%{(2, 1^2)), k^k+l 

4.2.1 - If U2 < Yi'il] '^2/(:+i(^2i 1^2) then with probability 1 /2 accept sample path, else return to step 1 

4.2.2 - If U2 > n*=i^ '^2k{^^^ ''^^^ reject sample path and return to step I 
5 - Discard or Retain Auxiliary Information 



Upon accepting a proposed sample path skeleton within the UEA (as simulated by Algorithm 6. 1 and Algorithm 
6.2 and so satisfying P1-P2), we need to be able to simulate the sample path at further inference times (Algorithm 
3.1 Step *) in order to satisfy P3. Any further simulation is conditional on information obtained constructing the 
sample path skeleton. In particular, our sample path belongs to D, (by Algorithm 6.1), the sample path minimum 
(or maximum) belongs to a particular layer (as a consequence of the mixture proposal in (52 1), we have simulated 
the sample path minimum (or maximum) (either Xr - ihsj or X^ = m,., by Algorithm 5.2) and skeletal points 
(X^i , . . . ,X^^) and finally we have simulated whether the sample path maximum (or minimum) is contained in the 
first (i - 1) layers or in the /,* layer (by evaluating the Radon-Nikodym derivative in (53 1 by means of (54i and 



p5\ ). In summary, we have four possible sets of conditional information for our sample path X, 

5 1 := ^Xs,X,,X e Di,msj e [(jcAy) — a,, (xAy) — fl,_i jjX^ = msj,X^^, . . . ,X^^, ihs^ e |^(jicV3'), (xWy) + fl,-ij} 
^2 := [x,,X,,X e Di,msj e [(xAy) - fl,,(xAy) - fl,-i],XT = ms,,,X^,, 

53 := ^Xs,X,,X e D,,msj e \(xWy) + fl,_i,(jcVy) + a,],^^ = ms,t,X^,, 

54 := i^Xs,X,,X e Di,msj e [(.xVy) + fl,_i,(jcVy) + fl,],^^ = mj,,,^^,. 



e [(.xVy) + fl,_i,(xVy) + ai]} 
e [(xAy)-fl,_i,(xAy)]} 
rhs,, e [(xAy) - a,, (xAy) - 



The difficulty in simulating further intermediate inference times conditional on the above is that the layer infor- 
mation pertaining to the sample path minimum and maximum induces a dependancy between the sub-interval in 
which we want to simulate an intermediate point, and all other sub-intervals. An additional complication arises 
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as we know precisely the minimum (or maximum) of the sample path, so the law we need to simulate further 
inference points from is that of a Bessel bridge conditioned to remain in a given interval. 



However, the minimum (or maximum) simulated in Algorithm 6.2 Step 2 is auxiliary sample path information 
(as in ( |46] l) and doesn't constitute part of the exact algorithm skeleton, so can be discarded. Furthermore, layer 
information regarding the sample path minimum and maximum is sufficient in determining a layer for the entire 
sample path. As such, reconsidering S i in light of the above (52, ^3, 54 can be similarly considered) we have. 

Si := \^X,,X,,X^^,...,X^^,ms,, € [(jcAy) - a^, (xAy) - fl,_i], m.,,, e [(xVy), (xVy) + a^^i]} 
Now, to remove the layer induced dependency between sub-intervals we can unbiasedly simulate for each sub- 



interval a layer for the sample path minimum and maximum within that sub-interval, as outlined in Section 7.3 
and Algorithm 7.3. Further intermediate inference points can then be simulated in a variety of ways as outlined in 
Sections |T2]|73] and |T6) Algorithm 7.5 and Algorithm 7.6. 



6.2 Localised Approach 

The Localised Approach is based on the layered Brownian bridge construction found in the Localised Exact Algo- 
rithm (LEA) originally proposed in IfTTlfTTl . The LEA is an alternative exact algorithm based on the mathematical 
framework of EA3 (see [TJ). However, we don't go into detail as to its construction as in the context of this paper 
it suffers from a number of computational weaknesses (in particular sigrrificant computation is required in order 
to satisfy P3) so is not well suited for the purposes of our paper. Instead, we briefly outline its construction and 
highlight which aspects of its construction present difficulties. 

The key notion in the Localised approach is that rather than proposing sample path skeletons from K ' (where 
the end point Xt '■- y ~ h is first simulated), the interval to be simulated ([0, T]) can be instead broken into a 



number of bounded segments (as in (18 1). Each segment is successively simulated by means of simulating the first 
hitting time r, of a Brownian motion proposal sample path (as outlined in 1 11 1) of some user specified boundary 
symmetric around its start point (for instance, ifX^-x with boundary 9 then t := inf{5 : Xs i {x - 6,x + 9]]), and 
simulating and accepting a sample path skeleton conditional on the simulated boundary (with a suitable modifica- 
tion of the acceptance probability to account for the modified proposal measure). 

The benefit of the Localised approach is that simulating the first hitting time of a boundary acts as a layer for 
the bounded segment (i.e. Vm e [0,t],X„(w) e [x- 9,x + 6}) and so <P{Xq j) is conditionally bounded (as per R3) 
and furthermore a bound can be found for A(Xr) in the acceptance probability. However, as with the UJEA and 
AUJEA this approach to simulating sample path skeletons can result in simulating skeletons for intervals exceed- 
ing that required (which is computationally wasteful), further complicated by the need to specify the boundary 0. 
Furthermore (as discussed in |20|), this methodology can't be used to simulate conditioned diffusion and jump dif- 
fusion sample path skeletons (sample paths conditioned to hit some specified end point). Finally, unlike the Bessel 
approach, the minimum or maximum that is simulated forms part of the skeleton and so can not be discarded. As 



such, the demarginalisation strategy taken in Section 6.1 in order to extend the UFA with the Bessel approach for 



simulating layered Brownian bridges to satisfy P3 can't be conducted, so the skeleton can't be augmented. 



7 Adaptive Layered Brownian Bridge Constructions 



In Section p^ we proposed the AUEA (Algorithm 3.2) as an alternative to the UEA (Algorithm 3.1). In this sec- 
tion we outline how to simulate finite dimensional skeletons of layered Brownian bridges for use within the AUEA 
(and by extension the BJEA (Algorithm 4.1) and AUJEA (Algorithm 4.3)). In particular, we present new results 
for simulating initial (intersection) layer information (Algorithm 3.2 Step 2 - Section [TTT) , intermediate points 
conditional on the layer information (Algorithm 3.2 Step 3.1.2 - Section 7.2 1 and finally, new layer information 
for each sub-interval created by the intermediate point (Algorithm 3.2 Step 3.1.4 - Section[73]l. 



We use the results we present in Sections 7.1 7.3 to outline two novel layered Brownian bridge constructions 



which can used within the AUEA. In particular, in Sections 7.5 and 7.6 we introduce the Intersection Layer Ap- 



proach and the Hybrid Approach respectively, both of which ensure the AUEA satisfies P1-P3. 
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7.1 Simulating an Initial Intersection Layer 



Upon simulating a proposal Brownian bridge layer as per Algorithm 6.1 in Section 6.1 we know that our entire 
Brownian bridge sample path is contained within the i* layer, but is not contained within the (i - 1)* layer 
Simulating sample path intermediary points is complicated by this conditional information (and as discussed in 
Section|6j it is not possible to simulate intermediary points directly). The novel approach we take in this paper is 
to simulate further layer information regarding the minimum and maximum of the proposed sample path (which 
together provide a sample path layer). To achieve this recall (with reference to Section |3] and ( 49|50|5l" l) that. 



having simulated a layer for our proposal Brownian bridge sample path as per Algorithm 6. 1 , we know the sample 
path belongs to the set of sample paths D,. We can then simply decompose the set D, into a disjoint union of 
possible sets our sample path belongs and simulate (unbiasedly) which our sample path belongs to, 

A = -L, u f/, = (L, n f/,) w (t/f n L,) w [if n f/,), (56) 



This decomposition can be interpreted as the sample path attains the layer at both its minimum and maximum 
(Dj i) or its minimum (Z),,2) or its maximum (Di j,). We can simulate unbiasedly which set our sample path belongs 
to (which we call an Intersection Layer) by application of the following novel results and Algorithm 7.1. 

Theorem 7.1 (Initial Intersection Layer). The probability a Brownian bridge sample path is in D, i, given it is in 
D„ can be represented as follows (denoting fi :— (xAy) — Qj, fl :— (xAy) — a,_i, vi :— (xV3') + fl,_i, vl :— (x\/y) + aj, 



PD,, :=p(a,i|A,W, = x,H^,=3;) 

I3{s,t,x,yji,{'\,vi,v'\) 



P (s, f , X, y, ei, vi, v'\) + p (s, f, x,y, ti, £t (xWy), vl) + /3 (s, t, x, y, {x^y), vi, uT) 



(57) 



Proof. Follows by Bayes rule. Theorem 5.6 and the decomposition of D, in (56 1. □ 



Corollary 7.1. Events of probability po^ , can be simulated by retrospective inversion sampling ( as per Corollaries 



5.1 5.2 and Algorithm 5.1), noting that po^ ^ is a function of/3 probabilities, using the following sequence. 



Sl(s,t,x,y,n,{tvi,v'[) 



' Sf^^is, t, X, y, n, n, y^) + Sf^, (s, t, X, y, U, n,(xWy), vi) + Sf^^ (s, t, x, y, n,(xAy),vi, uT) ' 
Remark 7.1. By symmetry we have that P(A,2 I A, W, = x, W, =y) = P(A3 I A, Ws = x, W, =y). 



Algorithm 7.1 - Simulation of an Initial Brownian Bridge Intersection Layer 

1 - Simulate layer I — i as per Algorithm 6.1, simulate u ~ U[0, 1] and set j — 
2- Whileue{s^f''\S^j^l'),j^j+l 

3 - Ifu< S"^''^^ then set A = A,i 



4 - If u> then with probability 0.5 set A — A,2 ^1^^ set A — A3 



7.2 Simulating Intersection Layer Intermediate Points 

Having simulated an intersection layer we require a sampling scheme for simulating the conditional Brownian 
bridge at some intermediate time q € is, t). As shown in IH, the density at q can be written as follows, 

n(w) P(W, = w\w,, W„m,, e [4^, m.„ € [u^,,, v]^,]) (x p[s, q, t, x, w,y, v^, v^) ■ n{h„, ai). (58) 

A computationally efficient (though inexact) method of simulating from 7r(w) was outlined in [8 1 based on in- 
version sampling and numerical methods. Here we show that it is possible to extend |^, and simulate from this 
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density exactly by means of composition sampling (see f32i) and rejection sampling. We will begin by considering 
the upper convergent bounding sequence of p (ik e S>o we have p < S^i.^^ and Mmk^oo 'S2i:+i ~ P^- Considering 
the decomposition of ( [43] l into its elementary form in terms of <; and (f (see ( [33] l and (|34]| respectively) note that it 
is composed of K - 64(fe + 1)^ of these elementary terms. 

Recalling that <;';^,{j\x,y) := g;"(j\x,y) + s-;f'""(;; -.x, -y) and (fil-"(j-x,y) := x,y) + ^;J'"''(;; -a;, -y) it 

can be shown that each of these terms has the structural form exp {a, + biw}. As such, we can find a bound for our 
target density ( [58] l as follows (c, = {0, 1 ) denotes whether the density contribution is positive or negative), 

K 

n(w) oc p(w) ■ N cri) < S^.^^ ■ N al) = ^ [(-1)'' ■ exp {a, + biw) ■ 1 (w e [tt, Vi\) ■ N (rl)] 



= 2 f(- 1)'' ■ {«'■ + A'"'^'- + bicrl/2] ■ 1 (w e [^i, y,]) ■ A? + bicrl, crl)] 



1=1 

=: ^ o)/ ■ A? (a^/, (rf, ) ■ 1 (w e [^,-, u,]) . (59) 

!=I 

Here we have a mixture of positively and negatively weighted truncated Normal densities (with common variance). 
Although each truncated Normal in the mixture is unique (due to the truncation points), a significant number have 
common mean. For instance, if 2^ + 1 = 1 we have the Normal arising from gf^ifij; x, w) and g^'f {j; x, w) have 



the same mean. We exploit this by partitioning the interval that provides support for the target density ( 58 i into 
sections corresponding to the truncation points. As a consequence, the resulting mixture density has a number of 
positive and negative elements which cancel each other out (i.e. they can be netted from one another). Defining 
cu^ := {co V 0) we can find an upper bound by solely considering the mixture formed from the positive weights. 

K 

n{w) <Y,^i-N {fii, crl) ■ 1 (w e [£,, v^]) [l (w € [4,, + 1 (w e [<„ u^j) + 1 (w e [v[„ v],])] (60) 
1=1 

i 

^:S%\,-N crl). (61) 

By application of composition sampling (see |'32l) we can simulate from -Sj^^j ■ n(^iIk, o-^) by first choosing one 
of the truncated Normal densities partitioned on the interval [L, Y] with probability proportional to, 

uolj ■ [o (t \fi,, crl^ - (5 (l I^,., tr^, )] . (62) 

As w ~ 2j^j ■ (^fi„, (tI) and we require w ~ p ■ N (//„,, cr^) we accept this draw with probability, 

^ ^ p(w)-N{w\ti,,.,crl) ^ p(w) 

Events of probability P can be simulated by retrospective inversion sampling (as per Corollaries |5.1[ |5.2| and 
Algorithm 5.1), noting that P is a function of p. The complete rejection sampler is presented in Algorithm 7.2. 

Algorithm 7.2 - Simulation of Intersection Layer Intermediate Points 

1 - Simulate u ~ U[0, 1] and set j — I 

2 - Simulate w ~ 5^'^ ■ n(/Iw, o"w) for some k G S>o 



2 - While u e 



^5^/w) 5^,,,(w)N 



5^'^w)' 5f^w) j 
w) 



3 - If u < — -f then accept else reject 
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Simulating intermediate points as per Algorithm 7.2 is (typically) highly efficient as 5^'^ ■ N(ji„, cr^) (typically) 
forms a tight bound to jt(w) (as noted in |8|). If this is not the case, then sampling from the bounding density 
with k > 1 isn't usually effective as ^^^^j is only formed by the positive netted components of ^j^^j. We present 



an alternative (exact) sampling scheme in Section 7.6 which can be used in this case based on the Bessel Approach. 



Alternatively, it can be shown that as our target density n{w) < ^j^^j ■ A^(;U„,, crf^.) is constrained to the inter- 
val [^j,, {],], then by sampling the density S'j^^j ■ N(jj.h., cr^^,) on a sufficiently tight mesh over the interval [f^^, 
and applying a convexity argument (noting the bounding density is a mixture of Normal densities), we can find 
a uniform bound for S'jj^i ■ Nifin, cr^^,) for any k (and so a uniform bound for 7t{w)). This uniform bound can be 
used to construct a rejection sampler in the case where 5^'^ ■ Nifi^, crlj) doesn't form a tight bound to n{w) and we 
require a bounding density with A: > 1. Furthermore, we conjecture that n{w) is at most bimodal (with at most one 
mode in (^J, + yj,)/2] and at most one in [{(1, + v\.,)l2, Wj,]). Together these lead to highly efficient search 
methods for sampling exactly from ji(w). 



7.3 Dissecting an Intersection Layer 

Upon simulating intermediary points of a Brownian bridge sample path conditional on an intersection layer (for 
instance in Section 7.2 1, simulating further intermediary points in a sub-interval between any two existing consec- 
utive points is more complicated as there is a dependancy between all sub-intervals (induced by the intersection 
layer). To simplify this problem we can dissect an intersection layer into separate intersection layers for each pair 
of consecutive points by considering all possible dissections and unbiasedly simulating which one of these occurs. 



To guide intuition we first consider the case where we have a single intermediary point {Wq - w) within an 

IT IT 

existing intersection layer {Ws - x,Wt - y,fns,t e [{,,,£],i],ms,t 6 [i/J,, l»^,]) and we want to simulate separate 
intersection layers for the intervals [s,q] and [q, t] conditional on the known intersection layer and the simulated 
point. We begin by noting that the simulated point provides further detail on the interval in which the minimum 
and maximum lies. In particular, if w e ,] we have that mj , e [{\ w] and similarly if w e [y^ ,] then 

T T* T I * I 

we have that OTj_, e [w, v^,]. As such we denote :— ((\i,f\ w), v\, :- (v].,\/ w) and we now have, 

D.KWg = w) = {m„ e [4„ f;,]} n {m,,, e [vi;, vl,]}. (64) 

The attainment of a particular layer in the interval [s, t] by either the minimum or the maximum implies that the 
same layer is attained by the sample path in at least one of the sub-intervals [s, q] or [q, t]. As such, in our case 
there are 9 possible (disjoint) bisections (which we denote as Bi-Bg where B := Dt\{Wq - w) - iiJ^?^jB,) as 
illustrated in Figure |7]l. For instance, our sample path may lie in B(„ which more formally has the form, 

B, := {{m.,, e [<; (xAw)]} f] {m.,, e [vi;, vl]}} [j {{m,,, e f] {m,, e [(wWy), vi;]}} . (65) 

This notion can be extended to the case where we have multiple intermediate points ('W := [Wg^ - w-i,. . . , Wg^ - 
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Figure|7]l: Illustration of 9 possible (disjoint) bisections 
w„)), and want to dissect the interval into separate intersection layers. In particular, we are dissecting a single 
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intersection layer into (n + 1) intersection layers, each with a layer for the minimum and maximum in their own 
sub-interval. As the sample path minimum and maximum must exist in one of the intersection layers there are 
b := (2*"+'' - 1)^ possible dissections B", . . . ,B'y We can simulate unbiasedly which of these dissections our 
sample path lies in by apphcation of the following novel results and Algorithm 7.3. 

Theorem 7.2 (Intersection Layer Dissection). The probability a Brownian bridge sample path is in B" conditional 
on B and 'W is as follows ( denoting £,{i) and U{i) as the lower and upper layer sets for B"), 



Proof. Follows directly by Bayes rule, Theorem 5.5 and Theorem 5.6 



(66) 



□ 



Remark 7.2 (Intersection Layer Bisection). In the particular case where we have a single intermediary point 
then the probability a Brownian bridge sample path is in Bj (conditional on B and Wq — w) reduces to that in [i8\l 
(denoting {^^fq, il'!q v\''q, v\''q and t^',, t^'', v^',, as the bounds for B,- in the interval [s, q\ and [q, /] respectively). 



Pb, 



P [s, q, X, w, 4-;, 4i, i^!q, v%) ■ p [q, t, w,y, t^, ^j, v^l,, v^^ 



(67) 



p[s,q,t,x.w,y,l\^,,t%v%v\) 

Corolla ry 7.2. Events of probability ps, can be simulated by retrospective inversion sampling (as per Corollaries 
and Algorithm 5.1), noting that it is a function of ^"^/3f^,'^J^^'\Q,'W) and^"^pl^,f^y'^^''^\Q,'W) probabilities. 



5.1 



5.2 



using the following sequence. 



gB(n.i) ._ 



5f (s, t, X, y, Q, nV, £(i), 1/(/)) 



(68) 



S';f,{s,t,x,y,Q,^Vji,jl,vi,A,y 
Remark 7.3. Unbiased simulation of the dissection the sample path lies in can be conducted by inversion sampling 



and a Cauchy sequence representation of the cdf of B {69 \. In particular, by application of Corollary 5.2 
sample path lies in B" if for some k > and u ~ U[Q, 1] we have u e (Cj^*"']' '\ Cj^"'"''). 

;=1 



Algorithm 7.3 - Dissecting an Intersection Layer 



1 - Simulate u ~ t/[0, 1] and set j — 1 and k — 

2 - While u 6 (Cff Cf^'^f ), k^k+l 

3-Ifu< C^^'^^ set dissection layer B - Bj else set j - j + I and return to step 2 



7.4 Refining an Intersection Layer 

Suppose we have already simulated layers for the maximum and minimum of our proposal Brownian bridge 
sample path (m^,, e [^^,,4,] and m,,, e [v[^, l»J,]), but we require more refined layer information (i.e. we want a 
set of narrower layers \{\, - < - or \v\, - v\t\ - I^^L ~ '^sfl)- This can be achieved by noting that the 
sample path falls in one of the following 4 possible (disjoint) intersection layer refinements (where R :- Wi'^j/?,) 

T IT T IT 

and simply unbiasedly simulating which one our sample path lies in (denoting , 6 [^J s /] "si ^ [^^j p J)' 

R, = {m,,, e [4„ n {m,,, e [v^, vl]} , R, = {m,,, € [<„ n {m,, € [vi„ vl]} . 

In a similar fashion to Section |7.3| we can simulate unbiasedly which of the intersection layer refinements our 
sample path lies in by application of the following established results and Algorithm 7.4. 
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Theorem 7.3 (Intersection Layer Refinement |I8] Section 5.3]). The probability a Brownian bridge sample path is 
in R, conditional on R is as follows, 

Corollary 7.3 (|8, Section 5.3]). Events of probability p^. can be simulated by retrospective inversion sampling 
(as per Corollaries \5.1\ \5^ and A I go rithm 5.1), noting that it is a function of 13 probabilities, using the sequence. 



' ' Sl,{sj,x,yj[,,el,vl,v],)l' 



(71) 



Remark 7.4. By application of Corollary 5.2 unbiased simulation of the refinement the sample path lies in can 
be conducted by inversion sampling and a Cauchy sequence representation of the cdfofR ( 72 1. In particular, our 
sample path lies in Rj if for some k > and u ~ U[Q, 1] we have u € (Cj^^j'', C*^^'''). 



C 



(72) 



Algorithm 7.4 - Refining an Intersection Layer 



1 - Simulate u ~ U[0, 1] and set j — I and k — 

2 - While u e (cf/\ C'^lfi), k = k+l 

3-Ifu< 'Sjj''' set layer R = Rj else set j — j + 1 and return to step 2 



7.5 Intersection Layer Approach 

The Intersection Layer Approach for constructing a layered Brownian bridge is a direct application of Sections 



7.1 7.3 In particular, we simulate initial intersection layer information for the sample path (Algorithm 3.2 Step 
2) by application of Algorithm 7.1. Now, in Algorithm 3.2 Step 3 we iteratively simulate skeletal (intermediate) 
points, then new intersection layer information conditional on the skeletal points. This can be achieved directly 
by application of Algorithm 7.2 and Algorithm 7.3 respectively. 

We present the iterative AUEA Step 3 in Algorithm 7.5 which can be additionally used to conduct Step *. S 
denotes the set containing all intersection layer information. The set is composed of (« - 1) elements correspond- 
ing to the intervals between n existing time points. In particular, each element (Saj,) between two successive time 
points (a < b) contains information regarding the sample path at the end points and an upper and lower bound for 
both the minimum and maximum of the sample path in that interval [Saj, '■— [a, b, Xa,Xb, ^, ^, v^^ ^, vl ^}). 

Algorithm 7.5 - Layered Brownian Bridge Simulation (Intersection Layer Approach) 

1 - For each intermediate point required (q) 

1.1- Select the appropriate existing intersection layer Sa,b from S such that q e (a,b) 

1.2- Algorithm 7.2 

1.3 - Bisect interval as per Algorithm 7.3 to find new intersection layers Sa.i] and Sq^h 
1.4- SetS = Su[S,,q,Sq,b}\Sa,h 
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7.6 Hybrid Approach 



The Hybrid Approach for constructing a layered Brownian bridge is simply a modification of Algorithm 7.5 Step 
2 (for reasons discussed in Section [7!2] l. An alternative is to simulate a single intermediate point using the Bessel 
approach outlined in Algorithm 6.2. However, a modification has to be made to the acceptance probability as 
the intersection layer provides more precise information regarding the interval in which both the minimum and 
maximum is contained. In particular, if we have simulated layer D, i then with probability 1/2 we propose the 
auxiliary minimum (else maximum) in the layer and then only accept the proposal sample path if the sample 
path maximum (else minimum) is contained between the (t - 1)* and i* layer. However, if we have simulated 
layer D,,2 (else D,,3) then we propose the auxiliary minimum (else maximum) in the i* layer and then only accept 
the proposal sample path if the sample path maximum (else minimum) is contained within the (i - 1)* layer. 

8 6-Strong Simulation of (Jump) Diffusions 

In this section we outline a novel approach for simulating upper and lower bounding processes which almost surely 
constrain (jump) diffusion sample paths to any specified tolerance. We do this by means of a significant extension 
to the e-Strong Simulation algorithm proposed in fS]. We consider applications of our extension, which are reliant 



on the new methodology presented elsewhere in this paper, in Sections 8.1 and 8.2 along with simulation results. 



As originally proposed in |8| and presented in Algorithm 8.1, e-Strong Simulation is an algorithm which sim- 
ulates upper and lower convergent bounding processes {X^ and X^) which enfold almost surely Brownian motion 
sample paths over some finite interval [0, T}. In particular, we have V m e [0, T] and some counter n, 

Xiin + 1) < Xt(n) < X„ < Xl(n + 1) < 4(«), (73) 

It should be noted that Algorithm 8.1 is an extension of that originally proposed in [SJ. In particular, in contrast to 
lID, we can now simulate an initial intersection layer (Algorithm 8. 1 Step 2) and simulate the intermediary points 
exactly (Algorithm 8.1 Step 3.1). 

Algorithm 8.1 - e-Strong Simulation of Brownian Motion sample paths (n bisections) 

1 - Simulate Xj :-y ~ N(0, T) and set / = 1 

2 - Simulate initial intersection layer S :— Sq t — \0, T , Xq, Xj , j., ij^ j., j, as per Algorithm 7.1 

3 - While i < n 

3.1 - For each of the 2'^^ intersection layers in S (denoted S{ , :— ^s-',t-',x{,xj,il'\,ij],vj\,v'fi'^) 

3.1.1 - Simulate Xg where q :- (s-i + f-')/2 conditional on Si., as per Algorithm 7.2 

3.1.2 - Bisect si f into Ss^q andS^^i as per Algorithm 7.3 

3.1.3- For Sj^g and S^^^,, while ^ej*'^ -ii*'K^> yjiti - si) I A or - vi*'K | > yjit' - si) I A then refine 
intersection layer as per Algorithm 7.4 

3.2 - SetS:^ U5=I [S^li U and i^i+l 



The intersection layer information can be used to find the bounding processes in ( 73 1, 

2"-' 2"-' 



/■=! (=1 

Furthermore we can find upper and lower bounds for the path integral, 

2"-' 2"- 



X\n) Yj <r ■ (f - ^'); ^H«) := Yj ^'i ' " "'^ ^^^^ 



1=1 !=1 
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It was shown in |l8l (and summarised below) that the dominating processes converge almost surely in the supre- 
mum and Li norms, where the rate of convergence in Li is of the order (9(n '^^). 



T I 

Theorem 8.1 (Convergence [8 , Prop. 3.1]). Considering X„(n) and X,i(n) as defined in ( 74 ), then convergence in 
the supremum norm holds in the limit ai « — » oo. 



w.p. 1: Hm sup |Xj(n) - Xi(n)\ 0. 
Theorem 8.2 (Li distance (H Prop. 3.2]). Considering the distance, 

\x^-x\ = f\xl-xil d« 

Jo 



(76) 



(77) 



then. 



2"/^ xE[|xT -X^IJ =0(1) (78) 
Now, considering the e-Strong Simulation of Jump Diffusions, note that upon simulating a jump diffusion sample 



path skeleton (as per the AUJEA), it has a form (see (22 1 and (28 1) that can be used in Algorithm 8.1. As such. 
Algorithm 8.1 can be extended to jump diffusions (Algorithm 8.2), and Theorems|8.1|and|8.2|still hold. 



Algorithm 8.2 - e-Strong Simulation of Jump Diffusion sample paths (n bisections) 

1 - Simulate jump diffusion skeleton as per Algorithm 3.2 to obtain initial intersection layer 

2 - Simulate further intersection layers as required (n bisections) as per Algorithm 8.1. 



As far as we are aware there does not exist any other (exact or inexact) method for the e-strong simulation of 
jump diffusions. The class of jump diffusions this methodology can be applied to is broad (the conditions outlined 
in Section |2] are sufficient) and motivate a number of avenues for future research and application. In particular, 
non-trivial characteristics of the diffusion path can be simulated (for instance extrema, hitting times, integrals) 
and can be applied to areas such as option pricing and the simulation of stochastic volatility models (which are 
currently being explored in related work). The precise implementation of Algorithm 8.1 and Algorithm 8.2 can 
be tailored to the specific application. 



8.1 Application 1 

In this section we consider the e-strong simulation of jump diffusion sample paths which can be represented as 
solutions to the following SDE, 



dX, = -X,. dt + dW, + dr 



Xo = 2, f€[0,5], A(X,.) ^ sm(X,X v(X, ) ~ f,(X,_) ^ N(-X,J2,l) (79) 



In this case, as the jump intensity can be bounded, we can simulate sample path skeletons satisfying P1-P3 us- 
ing an AUEA incorporated within the BJEA. In particular, for this example we have that (piX,) :- (Xj - l)/2, 
4>(Xr)\(Lx,Ux) G [-l/2,(max{(Lxf,(Uxf} - l)/2] and AiX,) < 1 A. Recall that in the BJEA the inter- 
val the sample path is to be simulated over ([0,5]), is broken into segments corresponding to proposed jump 
times C^i, . . .). As such, considering the simulation of a diffusion sample path in the interval [*I'i,*I'2] condi- 
tional on Xij<(i) then the proposed end point is simulated according BBM, ^4/(2) ~ h{X^^2y,X'jni-f,^'2 - ^\) 
exp {-X2(2)/2 - (X^(2) - X>p(i))' / [2(>I'2 - Ti)]}. 



In Figure 8.1 1 we consider a sample path skeleton simulated from the measure induced by (79 1 using the AUEA 



and BJEA. In Figure 8.1 2 we consider the e-strong simulation of the simulated sample path as per Algorithm 8.2 



with an increasing number of bisections. Finally, in Figure 8.1 3 we consider the e-strong simulation of a number 
of other sample paths. 
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Time 



Figure 8.1 1: Sample path skeleton generated from the measure induced by (79i. Vertical blue lines are the start 
and end point of the skeleton, vertical red lines are accepted jump times, horizontal black lines are almost sure 
bounds constraining the sample path and horizontal dotted red lines concern the layer attained by the sample path 
minimum and maximum in any given sub-interval. 



8.2 Application 2 

In this section we consider the e-strong simulation of jump diffusion sample paths which can be represented as 
solutions to the following SDE, 

dX, = sin(X,Odf + dW, + djf'\ = 0,f e [0,2],A{X,.) = xl,v(X,.) ~ fy{X,.) = U[-X,.aO,-X,.wO] (80) 

In this case the jump intensity is unbounded so to simulate sample path skeletons satisfying P1-P3 we use the AU- 
JEA. In particular, for this example we have that (p(X,) := (sm^(X,) + cos(X,))/2 G [-1/2,5/8], A(X,)\ (Lx, Ux) < 
max {(Lx)^, (Ux)^^ and the end point simulated according to BBM, Xj :- y ~ h oc (exp |- cos(y) - 3'^/6}. 



In contrast with Section 8.1 in this section we instead consider the e-strong simulation of jump diffusions gen- 
erated from the measure induced by ( [79| to some specified tolerance e (such that X'\ - Xi < e). As such we 
appropriately modify Algorithm 8.1 and Algorithm 8.2 and instead bisect intervals in which the upper and lower 
bounds for the sample path are less sharp. 



In Figure 8.2 1(a) we consider a sample path skeleton generated from the measure induced by ([79]l and in Figure 



8.2 1(b) we consider the e-strong simulation of this skeleton to the tolerance e = 0.5. In Figure [8^ 2 we consider 



the e-strong simulation of a number of other sample paths to variety of tolerances. 



28 



1 2 3 

(a) n = 1 
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(c) n = 4 
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(b) n: 



(d) n = 8 



Figure 8. 1 2: e-strong simulation of jump diffusion sample path in Figure 8.1 1 . n denotes the number of simulated 
bisections as per Algorithm 8.1. 
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Figure 8.1 3: Sample path skeletons generated from the measure induced by (79 1 with « = 6 bisections 
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(a) Sample path skeleton simulated from the measure induced by 
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(b) e-strong simulation of (a) to tolerance e = 0.05 



Figure 8.1 1: Vertical blue lines are the start and end point of the skeleton, vertical red lines are accepted jump 
times, horizontal black lines are almost sure bounds constraining the sample path and horizontal dotted red lines 
concern the layer attained by the sample path minimum and maximum in any given sub-interval. 
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Figure 8.2 2: e-strong simulation of jump diffusion sample paths from the measure induced by ( 80 1 



31 



References 



[1] Y. Ai't-Sahalia. Closed-form likelihood expansions for multivariate diffusions. The Annals of Statistics, 
36:906-937, 2008. 

[2] T.W. Anderson. A modification of the sequential probabihty ratio test to reduce the sample size. Annals of 
Mathematical Statistics, 31(1):165-197, 1960. 

[3] S. Asmussen, P. Glynn, and J. Pitman. Discretization error in simulation of one-dimensional reflecting 
Brownian motion. Annals of Applied Probability, 5(4):875-896, 1995. 

[4] O.E. Barndorff-Nielsen and N. Shephard. Power and bi-power variation with stochastic volatility and jumps. 
Journal of Financial Econometrics, 2(1): 1-37, 2004. 

[5] V. Benes. Existence of optimal stochastic control laws. SIAM Journal on Control, 9(3):446^72, 1971. 

[6] A. Beskos, O. Papaspihopoulos, and G.O. Roberts. Retrospective exact simulation of difiiision sample paths 
with apphcations. Bernoulli, 12:1077-1098, 2006. 

[7] A. Beskos, O. Papaspihopoulos, and G.O. Roberts. A factorisation of diffusion measure and finite sample 
path constructions. Methodology and Computing in Applied Probability, 10:85-104, 2008. 

[8] A. Beskos, S. Peluchetti, and G.O. Roberts, e-strong simulation of the Brownian path. Bernoulli, In Press. 

[9] A. Beskos and G.O. Roberts. An exact simulation of diffusions. Annals of Applied Probability, 15(4):2422- 
2444, 2005. 

[10] F. Black and M. Scholes. The pricing of options and corporate habilities. Journal of Political Economy, 

81(3):637-654, 1973. 

[11] Z.A. Burq and O. Jones. Simulation of Brownian motion at first passage times. Mathematics and Computers 
in Simulation, 77:64-71, 2008. 

[12] B. CaseUa and G.O. Roberts. Exact simulation of jump-diffusion processes with Monte Carlo apphcations. 
Methodology and Computing in Applied Probability, 13(3):449^73, 2010. 

[13] N. Chen and Z. Huang. Localisation and exact simulation of Brownian motion driven stochastic differential 
equations. Mathematics of Operational Research, In Press. 

[14] L. Devroye. Non-Uniform Random Variate Generation. Springer, 1st edition, 1986. 

[15] B. Eraker, M. Johannes, and N. Poison. The impact of jumps in volatility and returns. The Journal of 
Finance, 58(3): 1269-1300, 2003. 

[16] P. Etore and M. Martinez. Exact simulation of one-dimensional stochastic differential equations involving 
the local time at zero of the unknown process. ArXiv e-prints 1102.2565, 201 1. 

[17] K. Giesecke and D. Smelov. Exact sampling of jump-diffusions. Operations Research, In Press. 

[18] A. Golightly and D.J. Wilkinson. Bayesian sequential inference for nonlinear multivariate diffusions. Statis- 
tics and Computing, 16(4):323-338, 2006. 

[19] A. Golightly and D.J. Wilkinson. Bayesian inference for nonlinear multivariate diffusion models observed 
with error. Computational Statistics &• Data Analysis, 52(3): 1674-1693, 2008. 

[20] F.B. Gongalves and G.O. Roberts. Exact simulation problems for jump-diffusions. Methodology and Com- 
puting in Applied Probability, In Press. 

[21] M. Johannes, N. Poison, and J. Stroud. Optimal filtering of jump-diffusions: Extracting latent states from 
asset prices. Review of Financial Studies, 22(7):2259-2299, 2009. 

[22] I. Karatzas and S. Shreve. Brownian Motion and Stochastic Calculus. Springer- Verlag, New York, 2nd 
edition, 1991. 



32 



[23] J.F.C. Kingman. Poisson Processes. Clarendon Press, 1st edition, 1992. 

[24] RE. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations. Springer, 4th edition, 
1992. 

[25] R.C. Merton. Theory of rational option pricing. Bell Journal of Economics and Management Science, 
4(1):141-183, 1973. 

[26] R.C. Merton. Option pricing when underlying stock returns are discontinuous. Journal of Financial Eco- 
nomics, 3(1): 125-144, 1976. 

[27] B. 0ksendal. Stochastic Differential Equations. Springer, 6th edition, 2007. 

[28] B. 0ksendal and A. Sulem. Applied Stochastic Control of Jump Diffusions. Springer, 2nd edition, 2004. 

[29] U. Picchini, A. Gaetano, and S. Ditlevsen. Stochastic differential mixed-effects models. Scandinavian 
Journal of Statistics, 37(l):67-90, 2009. 

[30] E. Platen and N. Bruti-Liberati. Numerical Solution of Stochastic Differential Equations with Jumps in 
Finance. Springer, 1st edition, 2010. 

[31] K. Potzelberger and L. Wang. Boimdary crossing probability for Brownian motion. Journal of Applied 
Probability, 38:152-164, 2001. 

[32] B.D.Ripley. Stochastic Simulation. JohnWUey, 1st edition, 1987. 

[33] CP. Robert and G. CaseUa. Monte Carlo Statistical Methods. Springer, 2nd edition, 2004. 

[34] G. Sermaidis, O. Papaspiliopoulos, G.O. Roberts, and P. Feamhead. Markov chain Monte Carlo for exact 
inference for diffusions. Scandinavian Journal of Statistics, In Press. 



33 



